ABSTRACT 


LIN, CHUNG- YI. Determination of the Fracture Parameters in a Stiffened Composite 
Panel. (Under the direction of Dr. F. G. Yuan) 

A modified J- integral, namely the equivalent domain integral, is derived for a 
three-dimensional anisotropic cracked solid to evaluate the stress intensity factor along 
the crack front using the finite element method. Based on the equivalent domain integral 
method with auxiliary fields, an interaction integral is also derived to extract the second 
fracture parameter, the F-stress, from the finite element results. The auxiliary fields are 
the two-dimensional plane strain solutions of monoclinic materials with the plane of 
symmetry at *3=0 under point loads applied at the crack tip. These solutions are 
expressed in a compact form based on the Stroh formalism. Both integrals can be 
implemented into a single numerical procedure to determine the distributions of stress 
intensity factor and F-sfiess components, T\\, Tn, and thus F33, along a three-dimensional 
crack front. 

The effects of plate thickness and crack length on the variation of the stress 
intensity factor and F-stresses through the thickness are investigated in detail for through- 
thickness center-cracked plates (isotropic and orthotropic) and orthotropic stiffened 
panels under pure mode-I loading conditions. For all the cases studied, T\\ remains 
negative. For plates with the same dimensions, a larger size of crack yields larger 
magnitude of the normalized stress intensity factor and normalized F-stresses. The results 
in orthotropic stiffened panels exhibit an opposite fiend in general. As expected, for the 



thicker panels, the fracture parameters evaluated through the thickness, except the region 
near the free surfaces, approach two-dimensional plane strain solutions. In summary, the 
numerical methods presented in this research demonstrate their high computational 
effectiveness and good numerical accuracy in extracting these fracture parameters from 
the finite element results in three-dimensional cracked solids. 
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1 Introduction 


The study of fracture mechanics emerged in the early twentieth century. Among a 
handful of researchers, Griffith's idea of “minimum potential energy” [1] provided a 
foundation for all later successful theoretical studies of fracture, especially for brittle 
materials. But it was not until after World War II that fracture mechanics developed as a 
discipline. Derived from Griffith's theorem, the concept of energy release rate, G, was 
first introduced by Irwin [2], and was in a form that is more useful for engineering 
applications. He defined the energy release rate, or the crack extension force tendency so 
that it can be determined from the stress and displacement fields in the vicinity of the 
crack tip rather than from considering an energy balance for the elastic solid as a whole, 
as Griffith suggested. Irwin also used the Westergaard stress function [3] to show that the 
stresses and displacements near the crack tip of an isotropic linear elastic material in 
mode-I plane stress could be described by a single parameter, K , that is related to the 
energy release rate [4], i.e., 

G = K 2 /E, (1.1) 

where E is the Y oung's modulus. For plane strain, E is replaced by E / (1 - v 2 ) . This crack 
tip characterizing parameter later became known as the stress intensity factor. 

Rice [5] later defined a path-independent ./-integral for two-dimensional crack 
problems in linear and nonlinear elastic materials. As shown in Figure 1.1, J is the line 
integral surrounding a two-dimensional crack tip and is defined as 

J = J (Wdx 2 - (Tyfij ds ) , i,j= 1,2 (1 .2) 

r 1 

where r is a curve surrounding the crack tip, W is the stress-work density, <jjj are 
components of the stress tensor, tij is the j- th directional component of the unit normal 
vector on the path r, and ds is an element of arc length along F It was shown that the J- 
integral is a more general type of energy release rate. For a linear elastic material, G = J . 
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Therefore, the stress intensity factor K can be readily obtained, according to Eq. (1.1) and 
the computational efficiency of the ./-integral, as 

K=41e. (1.3) 

The /-integral is effective for evaluating K in two-dimensional crack problems. 
For three-dimensional problems, however, it is difficult to distinguish K at different q 
locations, assuming the line integral is performed on the X 1 -X 2 plane. Thus an alternative 
procedure needs to be developed to determine the distribution of K through the thickness. 
Parks [6] employed the virtual crack extension method to determine J from elastic -plastic 
finite element solutions. The method is based on an energy comparison of two slightly 
different crack lengths and requires only one elastic-plastic finite element solution, 
because the altered crack configuration is obtained by changing nodal positions. The 
procedure is directly applicable to two-dimensional configurations but can be extended in 
a straightforward manner to obtain arc-length-weighted average values of J along three- 
dimensional crack fronts. The three-dimensional applications, however, have significant 
loss of accuracy in the near-tip region where the values of field quantities (stresses, 
strains, and displacements) are required to determine the point-wise energy release rate 
along the crack front. Based on the virtual crack extension method, deLorenzi [7,8] 
developed a finite element method that is more general to calculate the energy release 
rate in two-dimensional and three-dimensional fracture problems and could include the 
effects of body forces and traction loading on the crack faces. 

Another investigation was made by Li et al. [9], They proposed a formulation 
which would convert area integrals to volume integrals in order to calculate point-wise 
values of the energy release rate along a three-dimensional crack front. Shih, Nakamura 
and co-workers [10,11] then developed a finite element formulation to calculate the 
volume domain integral. About the same time, Nikishkov and Atluri [12,13] applied a 
somewhat different approach but a similar numerical procedure, and named the 
formulation “equivalent domain integral (EDI)” which would be used by subsequent 
researchers [14-16], All of those derivations involve the application of the divergence 



3 


theorem and the implementation of a spatial weighting function that is based upon the 
virtual crack extension method. 

With the EDI method, a point-wise value of J along a three-dimensional crack 
front can be calculated, and therefore the value of K along the crack front can be obtained 
from Eq. (1.3). Another advantage is that the EDI method transforms surface integrals in a 
three-dimensional problem into integrals over a volume, or a domain (hence the name of 
equivalent domain integral), without evaluation of the stress singularities directly on the 
crack front. 

The stress intensity factor alone is not enough to characterize the crack behavior 
in some cases. Other fracture parameters may be needed to describe the crack behavior 
more precisely. As Irwin [4] pointed out there is a mathematical expression for crack-tip 
stress distributions in linear isotropic solids, Williams [22] showed that the expression is 
in fact an infinite power series of r, where r is the distance from the crack tip. The power 
series, in a concise form, can be written as 

a v (r,0) = , ij =1,2 (1.4) 

n =- 1 

where A„ are unknown constants which depend on the geometry and loading conditions, 
and fij n) (0 ) are the known angular distributions. The mode-I stress intensity factor is 
included in the first term of Eq.(1.4) as 

<r„=l m-Tr-AV). 0-5) 

r ^° 42m 

in which the stresses are singular at r = 0 and A ] = K l /427r ■ The leading term of the 
series of Eq.(1.4) represents r singularity; the second term is a constant; the third and 
higher-order terms are proportional to r (l/2) ", «=1,2,3.... Larsson and Carlsson [23] first 
denoted this constant term as T, and later it became the so-called “elastic 7-stress”. 

In addition to the stress intensity factor, the elastic 7-stress provides another 
parameter to identify the severity of stress and displacement fields near a crack. Larsson 
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and Carlsson [23] showed that the 7-stress is necessary to modify the solution of the 
stress state in a small-scale yielding crack problem in plane strain condition. Rice [24] 
showed that 7 is in fact a constant stress parallel to the crack flank, and included it as a 
second crack tip parameter to characterize suitably small plane strain yield zones. Several 
methods have been used to practically determine the 7-stress [25], In addition to the 
methods mentioned in [25], recently other methods were also used, such as the boundary 
layer method and the displacement field method [26], as well as the stress difference 
method [27], Among those methods, the interaction integral method developed by 
Nakamura and Parks [28] demonstrated highly computational effectiveness since it is 
based on the EDI method and has the capability to compute the 7-stress not only in an 
isotropic material but also in an anisotropic material. 

Under the NASA Advanced Composite Technology Program, Langley Research 
Center (LaRC) has performed fracture toughness tests for various types of wing structure 
specimens made from stitched warp-knit fabric composites. Variations of in-plane 
geometry and crack length were evaluated from three kinds of specimen geometry [29]: 
compact tension (CT) specimen with the crack aspect ratios 0.46 < a/w < 0.69 ; center- 
cracked tension (CCT) specimen with 0.26 < 2a /w < 0.42; single-edge notched tension 
(SENT) with 0.25 < a/w < 0.34 . 

Methods based on the equivalent domain integral and Betti’s reciprocal theorem 
were developed by Y uan and Y ang [29] to extract the fracture parameters - critical stress 
intensity factor and 7-stress. With the limited experimental data, the results tend to show 
that the critical mode-I stress intensity factor provides a satisfactory characterization for 
engineering applications of fracture initiation in the composite of a given laminate 
thickness, provided the failure is fiber-dominated and the crack growth follows in a self- 
similar manner. In addition, the high constraint due to high tensile 7-stress may be 
expected to inhibit the crack extension in the same plane and promote the crack turning. 


Recently, LaRC performed a mode-I test on a five-stringer panel manufactured 
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using the stitched warp-knit composite material. The crack initially extended in a self- 
similar manner and then turned parallel to the stiffener direction when the crack 
approached stiffeners (see Figure 1.2). In this dissertation, the effects of the geometrical 
attributes on the fracture behavior of this panel are investigated by using three- 
dimensional finite element analysis and linear elastic fracture mechanics to analyze the 
composites. Due to the high computational efficiency, the equivalent domain integral 
method is used to calculate the through-thickness K\ stress intensity factor and the 
interaction integral method is adopted to compute the through-thickness 7-stress 
components. The algorithms of the equivalent domain integral and interaction integral are 
implemented into a single computer program, which reads a set of finite element 
solutions from a given mesh as the input to determine the distributions of the fracture 
parameters along the crack front. The composites are modeled as linear, anisotropic, and 
homogeneous materials. For the purpose of verification and comparison, a similarly 
cracked plate structure without stiffeners is also analyzed with the same composite 
material properties as well as an isotropic material. 

The derivation of the EDI method is reviewed in Chapter 2 by the approaches 
mostly found in [15]. The derivation of the auxiliary fields necessary in the interaction 
integral method for an anisotropic material is presented in Chapter 3. Chapter 4 shows the 
procedure to determine the stress intensity factor and components of the 7-stress from the 
values of the equivalent domain integral and interaction integral. The finite element 
models used in this research are described in Chapter 5; the associated results are 
presented in Chapter 6. Finally, the summary and discussion is presented and suggestions 
for future research are made in Chapter 7. 
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Figure 1.2 The cracked stiffened panel to be analyzed in the dissertation. ( Courtesy of NASA Langley Research Center) 
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2 Equivalent Domain Integral (EDI) 


The derivation will assume a traction-free crack in a linear elastic material, with 
the intention of determining the mode-I stress intensity factor K\ through the thickness. 


2.1 Mathematical Formulation 


Let us consider a small cylindrical volume with radius e encompassing a segment 
of crack front of length A such that both e and A approach zero, as shown in Figure 2.1. A 
local coordinate system is defined so that the axes x\ and X 2 are perpendicular to the crack 
front, while x\ and Xj are lying on the crack plane. The volume is enclosed by five areas, 
namely the outer surface A e , two end surfaces A eJ and A&, the top crack surface A £ch and 
the bottom crack surface A^. 


J 


The local /-integral over the outer surface A e of the tube is defined as [30] 
du , 


= lim— [ 

A ->° A J 
£->0 * 


(Wn x - (J, '-AA-n Ad A . 
dx, 


i,j = 1, 2, 3 


( 2 . 1 ) 


In Eq.(2.1), W is the stress-work density, defined as W = ja j de ij , where Cg are 

components of the stress tensor, and % are components of the strain tensor. u t are 
components of the displacement vector; is the j - th directional component of the unit 
normal vector on the surface A e . Since this research will be limited only to linear elastic 
materials, the stress-work density is simplified as W = (o ij e jj )/2 . Note that 
displacements, strains, stresses are expressed in the local crack front coordinate system. 


For the purpose of simplicity in later derivations, let 
d u. 


0 = l Vn,-a,,^n r 


( 2 . 2 ) 


Then Eq.(2.2 ) can be rewritten in terms of boundaries shown in Figure 2.1 to complete a 
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surface integral as 


J = lim— 

A_>0 A 
£->0 


J QdA + J QdA + J QdA 


+Ap- 


Ap. t +Ap. 


(2.3) 


The evaluation of surface integrals in Eq.(2.3) is tedious and could lead to errors because 
singular terms on the crack front are included for numerical integration. Therefore, a 
modified form of the surface integrals is desirable, and this modified form would be the 
equivalent domain integral. 


Now consider two tubular surfaces, A and A e , as shown in Figure 2.2. A is an 
arbitrary surface enclosing A e on which the ./-integral is calculated. Aj and A2 are end 
surfaces connecting A and A e . (A-A e ) ct and (A r A) C b denote the top and bottom crack 
surfaces between^ and A f , respectively. An enclosed volume (V-V e ) is surrounded by all 
of these surfaces, which are called collectively A 1 , defined as 

Aj, = A- A e + (A- A e ) ct + (A e - A) cb + A l +A 2 . (2.4 ) 


Based on the virtual crack extension theory, deLorenzi [8] proposed the concept 
of virtual node shift that forms the definition of a spatial weighting function, which is 
called ^-function by some researchers [12-16,30], We will adopt this name throughout 
this dissertation and use the symbol s to represent this spatial weighting function. 

According to the configuration shown in Figure 2.2, an arbitrary but continuous ■S'- 
function is defined between A and A e so that the function has the following properties: 

s(x l ,x 2 ,x i ) = 0 on A, A e j and A ^2, A j and A2; (2.5a) 

6 ’( x 1 , x 2 , x 3 ) = s(x 3 ) on A £ . (2.5b) 


In order to complete the surface integrals over^z, the first integral in Eq. (2.3 ) can 
be rewritten as an integral over the closed surface A% and an integral over the physical 
boundaries ( A - A e ) ct + (A e - A) cb . And with the use of Eq.(2.5), Eq. (2.3) becomes 
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j* QsdA + | QsdA + j* QdA + j 0. 

(A—A £ ) c/ +(A £ —A) cb A £ i+A £ 2 A a 


1 sdA 


( 2 . 6 ) 


In Eq.(2.6) /is the area under the ^-function curve on surface A e and is defined as 

f = js(x 3 )dx 3 . (2.7) 

A 

The .v-function is equal to zero on both end surfaces of A e i and A&; therefore, 



= 0 and Eq.(2.6) remains as 






(A-A £ ) ct +(A £ -A) cb 


+ j OsdA . 

A ect +A ecb 


( 2 . 8 ) 


In Eq.(2.8) the negative sign of the first integral, which is over an enclosed 
domain, comes from the opposite direction of the outer normal vector to the surface A e of 
the volume (V-V £ ) in comparison with the normal vector to the surface A e in Figure 2.1 . 
The other integrals in Eq.(2.8) are actually on the crack surfaces. Therefore, we may 
separate integrals in Eq.(2.8) into a “domain” integral and a “crack face” integral, 
denoted as 


' J j- [i/i domain ^ ) crack face 


where (J) domain = - J QsdA , 

and (J ) crack face — 


J QsdA + J QsdA = J C/ rack fllcc sdA . 


(A-A^x+iAe-A)^ A +A e 


(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 


By recalling Eq.(2.2), Eq.(2. 10) can be written as 


U) 


domain 




dll: ^ 


J 

Wsn x - 

V 1 ) 

n J 


- 

- 


dA. 


( 2 . 12 ) 


The use of divergence theorem on Eq.(2.12) gives the following result: 
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GOd 


-J 


r ().s 3 m,. 3 .v \ 

IF <t,,— -- — f/F 

3 x, 3 xj 3 x 


_ f <V 


( d*i 


■ <X. 


3 x, 


MF . 


(2.13) 


Since the analysis is limited to linear elastic materials, it can be shown that the second 
integral in Eq.(2.13) is equal to zero [13]. Thus Eq.(2.13) is simplified as 


V) 


domain 


-l 


( 


v-r. 


W 


ds 


3m, 3.v 


3 


3 x, ij 3 x, 3 x 

v 1 1 J J 


\lV . 


( 2 . 14 ) 


On the crack surfaces, the first and third directional components of the unit 
normal vector n are equal to zero («, = n . = 0 ), according to the local coordinate system. 

The second component of n has the properties of « 2 =+l on the bottom face and 
n 2 = -1 on the top face. Upon substituting these conditions into Eq.(2.2), we have 


f 


ficrack face Wfl\ 


3 M, 


3m, 

<7 n — n 2 + 0 22 t 
3x, 3Xj 


/?2 ” 1 ” 


3m 3 

3 x, 


( 2 . 15 ) 


Since o 22 =o }2 =0 on the crack surfaces, Eq.(2.15) is then reduced to 


(^crack face ^*12 


3 mj 

3 x, 


(2.16) 


For a traction-free crack surface, o 12 =0. Thus the value of Q in Eq.(2.16) is equal to 
zero, and all integrals in Eq.(2. 1 1) vanish. 

Therefore, for a traction-free crack in a linear elastic material, the equivalent 
domain integral for the determination of A), in terms of displacements, strains, and 
stresses, can be conveniently expressed as 


J, 


= -!f 


w 


ds 


■ cr„ 


3m, 3 s 


3 x, u 3 x, 3 x 


W. 


( 2 . 17 ) 


To make a computable expression of Eq.(2.17), some numerical implementation needs to 
be defined. 
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2.2 Numerical Implementation 


The 20-node isoparametric brick-shaped elements are frequently used in the three- 
dimensional finite element analysis of linear elastic crack problems. The typical finite 
element mesh around the crack front is a fan-type mesh, as shown in Figure 2.3. The 
shaded area indicates a domain over which the equivalent domain integral is calculated. 
All elements in and beyond this domain are 20-node elements. The wedge-shaped 
elements attached on the crack front, however, contain only 1 5 nodes for each element. 

The ./-integral is the sum of the domain integral contributed by each element in 
the designated domain, e.g., the shaded area in Figure 2.3. That is, 

n e 

GO— = £(0),, (2.18) 

i = 1 

where ( J) i is the volume integral over the i-th element, and n e is the number of elements 
enclosed in the domain. 


In finite element modeling, the displacements are expressed by shape functions 
and nodal displacements: 

20 

k= l ’ 2 ’ 3 ( 2 - 19 ) 

j = i 

where AT = N ) is the element shape function for a three-dimensional solid 
element, and £, //, g are the element’s local coordinates that range between ±1. (u k )j is 


the displacement component at the /-th node where j is the local node number within an 
element. Then for the volume integral of a single element, Eq.(2.17 ) can be written as 


;ni 


w 


ds_ 


-cr 


j* 


diij ds 
dr | dx k 


iet[j y&TjdC . 


( 2 . 20 ) 


For an element with 2x2x2 Gaussian integration points, Eq.(2.20) can be modified in the 
form of numerical integration as 
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w '=-7 


z z z 

XIX 

m = 1 «=1 />=1 


rjr 9^ /T / 

(E- U v OS 


\ 


9xj 


|m’„ w 


- ,w det[j] 


( 2 . 21 ) 


In Eq.(2.21), W is the stress-work density, u' T is the vector of displacement derivatives, 


a is the stress tensor, s’ is the derivatives of the .v-function. and det[J] denotes the 
determinant of the Jacobian matrix. w rn , w„, and w p are integration weights, and they all 
have the values of unity for 2x2x2 reduced integration [31], 


Eq.(2.21) is the equation to be used for computation; therefore, the numerical 
implementation of each item in this equation needs to be explicitly expressed, as shown 
in the following sections. Once all items in Eq.(2.21) can be readily calculated, the J- 
integral over the domain can be evaluated from Eq.(2.18). 


2.2.1 The Jacobian Matrix 


The Jacobian matrix is defined by 

dx l dx 2 dx 3 

W W W 

_ dx, dx 2 dx 3 
dr/ drj drj 
dx { dx 2 dx 3 

Ji dc 

Each component of the matrix, according to 

dx k _ y dNj f 
3^“"^ Uk y ’ 

dx k Y dNj 

drj j-f dr] 

dx k _ dNj . _ 

d£-j-fdC {x )j ’ 


( 2 . 22 ) 

finite element theory, is defined as 

£=1,2,3 (2.23a) 

£=1,2,3 (2.23b) 

£=1,2,3 (2.23c) 


where (x k ) , is the local coordinate component at the y-th node within an element. 
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2.2.2 The Stress Tensor and Stress-Work Density 


The stress tensor a of a linear elastic material is a 3x3 symmetric matrix shown as 

(2.24) 

The stress-work density of the linear elastic material is (o ij £ ij ')l2 , or 


o = 


0-11 

( 7,2 

*^13 

( 7,2 

< 7 22 

*^23 

*^13 

*-*23 

< 7 33 


^ 2 (p"ll^ll +<722^22 + <7 33 £ 33 ) + <7i 2 £]2 + t7 23 £ 2 3 +<^’ 13 ^ 13 - (2.25) 

Note that Gy and £y are the stress and strain components from the finite element solutions 
output on the integration points. 


2.2.3 The Derivatives of the s-Function 


s’ is the vector containing derivatives of the 5-function with respect to the local 
coordinate system and is expressed as 


s = 


ds ds 
dx l dx 2 


ds 

dx x 


(2.26) 


To evaluate Eq.(2.26), the 5 -function must be defined first. Since the 5 -function is 
arbitrary and satisfies Eq.(2.5), it can be conveniently defined by the sums of the element 
shape functions as 


20 

f (£,/7,o=5> 


J S J ■ 


(2.27) 


j = 1 


For the 20-node brick-shaped element shown in Figure 2.4, the 5 -function is 
completely defined by specifying 5 |1C " = 5 |14 ' =1 and zero on all other nodes in order to 
satisfy Eq.(2.5). This definition yields a quadratic 5 -function over a single element. With 
the definition, Eq.(2.7) also can be evaluated and hence /= (2/3)41 where zl is the length 
of the domain segment [30], 
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It is clear that the ^-function is a function of the element coordinate system (q„ j), 
£), rather than the crack front coordinate system (xi, xz, xj,). Thus s’ should be expressed 
in terms of (£ rj, £) before it can be evaluated. This can be done by the chain rule, as 


shown in the following equation: 


ds 


ds dE, ds dlj ds diE, 


ds 

3x, 


dq 3xj dr) 3x, dC, 3x, 


dE, 

ds 

> = < 

ds dE, ds dr ) ds dE, 

= j- 1 

ds 

— 

dx 2 


dq dx 2 dr) dx 2 3 E, dx 2 


drj 

ds 


ds 3 1 ds dr) ds 3 C 


ds 

dx 3 


dq dx 3 dr/ 3x 3 d£ 3x, 


m 


J 1 is the inverse Jacobian matrix containing the following components: 



dn_ 


dx. 

dx. 

dx , 


dr) 


3x, 

dx-, 

3x, 


dr) 

dc 

3x 3 

3x 3 

dx 3 


(2.28) 


(2.29) 


The derivatives of the ^-function with respect to the element coordinate system, i.e. — , 


— and - — in Eq.(2.28), can be evaluated in the same manner as Eq.(2.23). 
dr/ dq 


2.2.4 The Derivatives of the Displacements 


u'^ is the vector of displacement derivatives and can be expressed as 


The components in Eq.(2.30) are the derivatives of the displacements with respect to the 
local coordinate system. Similar to the derivatives of the ^-function, they should be 
evaluated in terms of the element coordinate system (£, £). With the use of the chain 

rule on Eq.(2. 19), each component of Eq.(2.30) can be obtained by 


du 2 du 2 du 3 
3x, 3x, 3x, 


(2.30) 
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d«* _ Y f dN J , dN J dri i dN J d C V. \ 
a*j ^3£3x, r) // 3x, 3£ 3x, | 


A: = 1,2,3 (2.31) 


y 'N -\ <■ 

In Eq.(2.31), — — . — — and — — are the components of the first row of the inverse 

3xj 3xj 3x, 


Jacobian matrix of Eq.(2.29); 


dNj 3 Nj 


and 


ac 


are the derivatives of the shape 


functions that can be readily computed. 
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3 Interaction Integral 


The interaction integral is necessary for extracting the elastic T-stress from an 
existing finite element solution. It is based upon the formulation of the EDI method as 
well as a superimposed auxiliary stress field. Kfouri [32] gave the auxiliary stress field 
that is the analytical solution corresponding to a point force applied to a crack tip and 
parallel to the crack surface under plane strain in isotropic solids. For a three-dimensional 
crack, the point force becomes a uniform line load that is applied along the crack front, as 
shown in Figure 3.1(a). This stress field is a function of r, the distance from the crack 
front, and 9, the angle from x\ axis toward X 2 axis; but it is independent of the crack front 
location * 3 . 


Nakamura and Parks [28] applied the auxiliary stress field with the interaction 
integral and successfully calculated the T-stress distribution along the three-dimensional 
crack front. The auxiliary stress field, however, is valid only for isotropic materials. For 
anisotropic materials, the corresponding auxiliary fields have been derived using Stroh 
formalism [34], 


Similar to Eq.(2.17) of the equivalent domain integral, the interaction integral for 
mode-I loading in a given domain may be expressed as 




a ds 

~ <7ij£,j d^ + 


f 




duf 


V 


a dllf 

“ + G ii 

OX, 3x, 


ds 


()x ; 


YIV 


i,j= 1,2,3 


(3.1) 


where cr‘ , , and u. are the components of the auxiliary stress, strain, and 

displacement fields, respectively. For the purpose of numerical integration of each 
individual element in a domain, Eq.(3.1) can be written similarly to Eq.(2.21 ) as 

" 2 2 2 “ 


(u =^\222 -* jk e* k ^ + iux° + »yy 


f 


\_m = 1 n = 1 p=\ | 


w ffl w„w^det[j] 


(3.2) 


In Eq.(3.2), o a and u a denote the stress tensor and displacement vector of the auxiliary 
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fields, respectively. e a jk are components of the auxiliary strain tensor. These entities are 
expressed in terms of the components of the associated tensor or vector as follows: 




_a 

a w 

_a 

<T U 

r 

os 2 

b 


o a 

= 

_a 

<T n 

_a 

®22 

_a 

*^23 




_a 

(Tj3 

_a 

*^23 

_a 

<t 33 



II 

to 


1 & 22^22 

(u a 

)' T = 

dll x 

du\ 

du 



1 


3 xj 

dx 


<■331,33 


■2 (e 


12^12 


'23^23 


13^13 


); 


(3.3) 

(3.4) 

(3.5) 


Quantities of Eq.(3.3) and Eq.(3.4) can be obtained by straightforward 
substitution of auxiliary stress and strain fields. Components in Eq.(3.5) can be computed 
similarly to Eq.(2.31) as 

M: 

dx x 


'x a 20 

_ y I 

dr. ZJ 


dNj dig | dN j dr] | dN j dC, 


M 


dg dx x dr] dx\ d£ dx\ 



k= 1,2,3 


(3.6) 


where u & k are components of the auxiliary displacement vector. All of the other items not 

associated with the auxiliary fields are calculated exactly in the same way as the 
equivalent domain integral is. 


Since the auxiliary strain and displacement fields are derived from the auxiliary 
stress field which is a function of r and 0 , both are functions of r and #as well. All terms 
in Eq.(3.2), however, should be evaluated with respect to the local coordinates (xi, X 2 , X 3 ). 
Therefore, the auxiliary field calculation must be done by converting the Cartesian 
coordinates of nodes or integration points to the polar coordinates before substituting 
them into the auxiliary field formulation. The computation of the auxiliary displacement 
field is straightforward because it needs only substitution of all 20 nodes’ coordinates 
within an element. The auxiliary stresses and strains will need the coordinates of the 8 
integration points. The location of these integration points with respect to the element 
coordinate system is illustrated in Figure 3.2. 
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Let us recall the Williams expansion of Eq.(1.4) which can be generalized to 
three-dimensional problems. It is assumed that the asymptotic expansion of the stress 
field near the crack front location under general loading conditions can be expressed as 


3cm) 


°tj = 


e>(0)+T ii+ o( i). 


«=1(I) 


i,j= 1,2,3 


(3.7) 


where k\, hi, and /cm are local stress intensity factors, f,) n) (9 ) are angular distributions for 


the crack-tip field given by the two-dimensional deformation of anisotropic elasticity 
theory, and o( 1) represents other higher order terms. T,y are the non-singular 7-stresses, 
which have three distinct components, namely 

W- 


TOT 

-mi u 1 13 

0 0 0 
T u 0 7 


33 


(3.8) 


Tn is obviously the stress component acting parallel to the crack plane [24] and 
can be determined by the interaction integral with an imposed uniform line load f as 
shown in Figure 3.1(a). Similarly 7 b can be determined by using a different set of 
auxiliary fields. Instead of the line load perpendicular to the crack front and the X2-X3 
plane, a constant force f$ in X 3 -direction and through the full length of crack front should 
be imposed. This configuration, as shown in Figure 3.1(b), will yield an auxiliary stress 
field necessary to extract 7 b. 7 33 is a combination of T\\ and 7 b and can be readily 
obtained after the other two 7-stresses are determined (see Chapter 4). 

In the following sections, the derivations of both types of auxiliary fields are 
presented in order to determine all of the 7-stress components. 


3.1 Auxiliary Fields for T n 

In this dissertation, we will be concerned with the composite plate structures, 
which usually have at least one plane of symmetry in materials. The convention of 
orientation for a composite plate is that the plate is on the X 1 -X 2 plane while the X 3 is the 
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direction out of plane [34], Since most composite plates have at least one symmetry plane 
at X3=0, we will limit the derivation under this restriction. This kind of material is called 
the monoclinic material with the plane of symmetry at X3=0, or simply the monoclinic 
material about X3=0. 


The generalized Hooke’s law states the stress-strain relation in contracted notation 
as 


o c =C £ c , (3.9) 

where a — [p) <J 2 cr 3 <r 4 <r 5 <x 6 ] — [cr u <J 22 o 33 <7 23 <7 13 <r 12 1 (3.10) 

and c° = \e x e 2 e 3 e 4 e 5 e 6 J = [e u e 22 e i3 2e 2i 2e u 2e u J . (3.11) 

C is a 6x6 matrix, and is called the stiffness matrix in which the components Cy are 

material properties. A monoclinic material about x 3 =0 has the following form of the 

stiffness matrix: 



c 

'-'12 

c 

13 

0 

0 

C 16 

Cj 2 

r 

^22 

r 

'-'23 

0 

0 

Qe 

r 

'-'13 

r 

'-'23 

r 

'-'33 

0 

0 

C 36 

0 

0 

0 

C44 

Qs 

0 

0 

0 

0 

Q 5 

Q 5 

0 

C,6 

r 

'-'26 

c 

y -'36 

0 

0 

C 66 


The inverse of the stress-strain relation defines the compliance matrix s, as 


(3-12) 


£ C =S<7 C , (3.13) 

where s is the inverse of C. Thus the compliance matrix of a monoclinic material about 
X3=0 has the form of 


s = C' 1 = 


°13 

0 

0 

*16 


J 23 

0 

0 

*26 


°33 

0 

0 

*36 


0 

0 

0 

*44 

*45 


0 

0 

0 

*45 

*55 


0 0 


J 36 

0 

0 

*66 


(3.14) 


As stated earlier, the auxiliary fields for T\\ are independent of X3. This implies it 
is under the condition of two-dimensional deformation for which £3=0. By ignoring the 
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equation for 03 in Eq.(3.9), the stress-strain relation of the monoclinic material can be 
written as 

[a, <t 2 <j 6 er 4 aj = C4, e 2 e 6 e 4 ej, (3.15) 

where C° is the reduced stiffness matrix, shown as 

r c c c 0 0 1 

Mi M 2 M 6 u u 

Cj 2 C 22 C 26 0 0 

C°= C 16 C 26 C 66 0 0 . (3.16) 

0 0 0 c 44 c 45 

0 0 0 c 45 c 55 

The inverse of Eq.(3. 16 ) gives the definition of the reduced compliance matrix s°, as 

^11 S \2 S \6 0 0 

s \ 2 s ' 22 s ' 6 0 0 

4 4 4 0 0 . (3.17) 

0 0 0 4 4 

000 4 4. 



The components of s° can be also obtained by solving for Oj in the third equation 
of Eq.(3.13) that will yield 

1 6 

0-3 =0-33= ^3 (3.18) 

5 33 ^7 

Substituting Eq.(3.18) into the other five equations of Eq.(3.13) will have 



S i3 S 3j 

S 33 


i,j- 1, 2, 4, 5, 6 (3.19) 


According to Stroh formalism for two-dimensional deformations of an anisotropic 
elastic body [35], the characteristic equations have the reduced compliance as 
coefficients: 

s uM —2s 16 /j + (2s u +s 66 )^i — 2s 26 [1 + s 22 = 0 (3.20a) 
s' 55 ju 2 - 2s' 45 jU + s 44 = 0 . (3.20b) 

The solutions to Eq.(3.20 ) are the eigenvalues of elastic constants, ju a (oc= 1, 2, 3), where 
fix, i±i, Ji x , and Jl 2 are roots of Eq.(3.20a), and JU 3 , Z7, are roots of Eq.(3.20b). u (/ are 
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complex numbers, and ji a are the conjugates of ju a . 

Under the uniform line load /j shown in Figure 3.1(a), the auxiliary stresses are 
inversely proportion to r, or <r a °c r -1 . In Stroh formalism, the real form solution for the 

displacement u a and the stress function (p l due to the point forces can be written as 

2u a =- — I + S(0) L _1 f , (3.21a) 

n 

2(f>* = L(0)L _1 f , (3.21b) 

where S and L are Bamette-Lothe tensors, S( 0) and L( 6) are their associate tensors, f is 
the vector of the line load per unit thickness, and I is the 3x3 identity matrix. These items 


are defined as follows: 

f = L/l 0 0] T ; (3.22) 

S(d) = — Re{\.(ln(cos0 + // a sin^))B T }; (3.23a) 

71 

r\ 

L(0) = -- Re{B(ln(cos6> + ju a sin0))B T }; (3.23b) 

z, z 2 0 

IT 1 = j,', z 2 e 0 (3.24) 

0 0 (m.s ' u )”' 


The definitions of the terms in Eq.(3.23) and Eq.(3.24) are given as follows. 

For the purpose of simplicity, let =cos6+// a sin6. In Eq.(3.23), ( ) implies a 
diagonal matrix, thus 

In g 1 0 0 

(ln(cos^ + ju a sin^)) = 0 ln^ 2 0 . (3.25) 

0 0 ln^ 3 

A and B are 3x3 complex matrices containing Stroh eigenvectors and are defined as 



25 


A = 


kj> i 
Mi 

0 


kiPi 

k 2 q 2 

0 


0 

0 

^3(^45 — ^44/ Mi ) 


and 


B = 


— k l ju l -k 2 /i 2 0 

k { k 2 0 
0 0 —k 3 


where p\,pi, q\, qi can be obtained from 

Pa — ‘hi Pa ~ ^16 Pa ^12 = ^uMa ~ $ 26 *^22 / Pa • (X ~ 1 , 2 

k\, k 2 , h are normalization factors satisfying the following relations: 

2 k x (q x — p x p x ) — 1 , 2k 2 (q 2 — MiPi ) — 1 » 2/c, (s 44 f p 2 — s 44 ) — 1 . 

The components of L" 1 in Eq.(3.24) are defined by the following relations: 

Mi +Mi = Ti + zj ; 

M\ Mi — T 2 t z 2* » 
e = Tl Z 2 “JW 

m — (■S'44^55 — ^45 ^45 ) 


(3.26a) 

(3.26b) 


(3-27) 

(3.28) 

(3.29a) 

(3.29b) 

(3.29c) 

(3.29d) 


Upon substitution of Eq.(3.22) through Eq.(3.24), the stress function of 
Eq.(3.21b) can be obtained as 

Wi 


f=- 


f/n 


7T 


Ret (j) 2 f , 

k 

(3.30) 

r , (kf p [ \ng x + k 2 ju 2 \ng 2 )- z 2 (kf p x In g x + k 2 p 2 In g 2 ) 

(3.31a) 

(k x p x In g \ + k\p 2 In g 2 ) + z 2 ( k x In g , + k 2 In g 2 ). 

(3.31b) 


To determine the auxiliary stress field from the stress function, let t r be the 
traction vector on a cylindrical surface of r = constant which can be obtained as [33] 

r r d6 ’ 


(3.32) 
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or t = 


^Re( 

w 


where t r =z l (kfjufQ. l + k: juHl 2 )- z 2 (k l 1 /U ] £l ] + ICfjSl , ) 

and t = -z, (kf ju 1 Q l + klfiSl 2 )+ z 2 (k ] 1 £l 1 + lcil 2 ). 

In Eq.(3.34) Q a is defined as 

O rm = Aw _ -sin6+// g cos6 
ad cos 9 + jl a sin 0 


a= 1, 2, 3 


(3.33) 

(3.34a) 

(3.34b) 

(3.35) 


Then the stresses in the cylindrical coordinate system are 
o rr = [cos 6 sin 6 0]t r ; 

o r 3 = [0 0 1 } 


° 66 °r$ 


o gi — 0 . 


(3.36a) 
(3.36b) 
(3.36c) 

(T r i is found to be equal to zero after substituting Eq.( 3.33 ) into Eq.(3.36b). Substitution 
of Eq.(3.33) into Eq.( 3.36a) also gives < 7 rr , which is a function of r and 9 , as 

= •^ 1 ‘ S " Reji,, cos 9 + t r sin d\ 
nr 1 


a 


(3.37) 


It should be noted that the auxiliary fields in the calculation of the interaction 
integral of Eq.(3.1) are all in the Cartesian coordinate system. Hence the auxiliary 
stresses <7,* must be obtained by applying coordinate transformation on Eq.(3.36). The 

results, which will be able to be implemented into the numerical computing procedure, 
are given as follows: 

./Xi cos 2 9 


= 


(J 2 2 


(7 U 


nr 


- Re^,. cos 9 + C 2 sin d}; 


f/u sin 2 o 


nr 


Reji,. cos 9 + C 2 sin#}; 


f\s' u sin 6 cos 6 


nr 


Re{i ; , cos 9 + t r> sin #}; 


(3.38a) 

(3.38b) 

(3.38c) 


<3 = = 0 • 


(3.38d) 
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(733 can be obtained by substitution of Eq.(3.38) into Eq.(3. 1 8), which yields 


o 


33 


= _^;,b i cos^ +% sin i e) Re ^ ^ sine} 
7irs„ ' 


(3.38e) 


Once the auxiliary stress field is available, the auxiliary strain field can be readily 
obtained from the inverse of the stress-strain relation of Eq.(3.13). 


The auxiliary displacement field is readily available from Eq.( 3.2 1 a) which, after 
a series of substitution of Eq.(3.22) through Eq.(3.24), will yield u\ = 0 and non-zero 
items of u x and u\ as 


\ U \ 1 — fl S U \ Z 1 1 1 n ,. , 2 Re + Z 2^12 1 

KJ 2;r (U\ \z x S 2 x +z 2 S 21 I 


(3.39) 


where Sn, 621 , S 12 , and S22 are components of S (6) and are defined as 


S u = -kl/i x p x \ng x -klfi 2 p 2 \ng 2 ; (3.40a) 

S 1 2 = k\ p x In g l + k\ p 2 In g 2 ; (3.40b) 

S'jj =-k\/l x q x \ ng 1 - k 2 ju 2 q 2 \ng 2 ; (3.40c) 

S 22 = k\q x In g x + k\q 2 In g 2 . (3.40d) 


Without loss of generality, the magnitude of the line load may be chosen as f\=\ 
as f is arbitrary. This assumption, as well as the information on material properties and 
nodal coordinates, will enable the calculation of the auxiliary fields necessary for the 
evaluation of the interaction integral, which in turn will determine T\\. 

The auxiliary stress field for T\\ in an isotropic material was given by Nakamura 
and Parks [28] as 

< 7 “ = -—cos 3 6 ; (3.41a) 

nr 

a\ 2 = cos^sin 2 0 ; 


(3.41b) 
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cr* 2 = —— cos 2 9 sin 9 ; 
nr 

f 

o\ 3 = vcos 6\ 
nr 


(3.41c) 


(3.4 Id) 


°T 3 =On = °- 


(3.4 le) 

It can be shown from strain-displacement relations that the auxiliary displacement field is 


(i-y ; k 

kE 


In r + 


sin 6 

2(l-v) 


= _iL±il^[(i _ 2v)0 -sin^cos 0\ 
2nE 


u\ = 0 . 


(3.42a) 

(3.42b) 

(3.42c) 


3.2 Auxiliary Fields for T n 


The approach to determine the auxiliary fields for Tn is similar to that of T\\, 
except the line load f is replaced by a constant force fi, in the X 3 -direction, as shown in 
Figure 3.1(b). This configuration will change the line load vector f in Eq.(3.22) as 
f = [0 0 / 3 ]\ (3.43) 


which will yield a different stress function as 
0 

0 > - 

k\ In g. 


f =_A Re J 

nm 


(3-44) 


By applying a similar derivation following Eq.(3.32) and Eq.(3.36), the auxiliary stresses 
in the polar coordinate system can be obtained as 

£T r 3 =-^RefeQ 3 }, (3.45a) 

nmr 

o ge =o rg =o d3 =o tr = 0. (3.45b) 


Note that Q 3 is obtained through the definition of Eq.(3.35). 


The transformation of stresses from r-6 plane to X 1 -X 2 plane gives the results that 
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a\ x , o\ 2 , and <7 22 all are equal to zero because o ge = o rr = o r0 = 0 . The transformation 
of < 7 , 3 , and a® is conducted by the following relation [36]: 


Wn 

1*^23 


COS# 

sin# 


-sin# | <7 r3 1 

cosO J \(7 ei \- 


(3.46) 


This operation yields the non-zero auxiliary stress components 

<7j a 3 = - /jC ° s6 Re& 3 Q 3 ], (3.47a) 

nmr 

crt. = /;Sm6 RefcQ 3 }, (3.47b) 

nmr 

Then <r 33 can be obtained by the use of Eq.(3. 1 8) that also shows cr 33 = 0 . This operation 
assumes a monoclinic material about X 3=0 by using its compliance components shown in 
Eq.(3.14). 


The auxiliary strain field is also readily obtained from the inverse of the stress- 
strain relation of Eq.(3.13). And the auxiliary displacement field is available from 
Eq.(3.21a) as well. But with a different f, both u\ and u\ will be equal to zero while only 
u\ is the non-zero displacement as 

u\ = — (in r + Re {in g 3 }). (3.48 ) 

2 nm 


The magnitude of the constant force may also be chosen as f,= 1 because f is 
arbitrary. Elence the auxiliary fields can be calculated, and therefore the Tn for a 
monoclinic material about X 3=0 will be determined. 


For an isotropic material, the auxiliary fields may be obtained by applying the 
material properties on the stiffness matrix of Eq.(3.12). A similar derivation will yield 


*^13 


f 3 cos 6 
2 nr 


(3.49a) 


<7 


a 

23 


/ 3 sin 6 
2 nr 


(3.49b) 
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<1 = 0-22 = 0*33 = <^r 2 = 0 i 

fA i+v), 


u, = -- 


nE 


-In r . 


(3.49c) 

(3.50) 
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Figure 3. 1 Auxiliary line load on a three-dimensional crack: (a) uniform forces f\ 

normal to crack front: (b) uniform forces f parallel to crack front. 



Figure 3.2 Locations of the integration points inside an element. 
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4 Stress Intensity Factor and J-Stress 


The formulations regarding the equivalent domain integral and the interaction 
integral are implemented into a FORTRAN computer program, which uses the finite 
element solutions computed from another ANSYS program as the input. Once the 
equivalent domain integral and the interaction integral are evaluated, the stress intensity 
factor and T-stresses can be determined. In the following sections, formulation to 
determine both parameters will be shown in terms of those integral quantities and 
appropriate material properties. These formulations can be easily implemented into the 
FORTRAN program as well. The full contents of the ANSYS and FORTRAN programs 
are provided in the Appendices. 

4.1 Stress Intensity Factor 

For an anisotropic cracked solid, the energy release rate G is related to the stress 
intensity factor through [29,33,37] 

G=| k T L'k, (4.1) 

where k 1 = [/c n k l k m ] are stress intensity factors and L 1 is the inverse of one of the 

Bamett-Lothe tensors as shown in Eq.(3.24). In elastic materials, the energy release rate 

G is equal to the value of /-integral. For a pure mode-I crack in an elastic material, k\\ = 

km = 0, and Eq.(4.1) will reduce to 
/ 

J ] (x 3 ) = h fk^x 3 ), (4.2) 

where Xj, is the crack front location defined in Section 2.1, and J\(xi) is value of the 
equivalent domain integral, as defined in Eq.(2.17), on this location. Therefore the local 
stress intensity factor k\(xG can be determined as 

2/i(* 3 ) 

V s' n e 


k l {x i ) = 


(4.3) 
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For an isotropic material, L " 1 has the form of 

1 0 0 
0 1 0 

0 0 — — 
l-v 

The Eq. (4.3) will become 
k 1 (x 3 ) = 

which is the plane strain condition as shown in Eq.( 1.3). 




(4.4) 


(4.5) 


4.2 T-Stress 


The T-stress, in general, includes three components, namely T\\, Tn and I33, as 
shown in Eq.(3.8). T\\ and 7 b should be determined from the evaluation of interaction 
integrals, while T 33 can be readily obtained after T\\ and Tn are computed. 

Let / I) and / 2) be the values of interaction integral of Eq.(3.1) when the 
superposed uniform load is f\ and fi,, respectively. For an anisotropic material, it can be 
shown that the following relation between 7-Stress and interaction integral exists [37]: 


'4 41141 

_*^15 ^55. 14 J 


7 qx 3 ) * 13 , 

— — 

f s a 
I (2 \x 3 ) 

A 


where £»(Xi, ) is the crack front extension strain at a given crack front location X3. / 1 '(xi) 
and / 2) (x 3) are the interaction integrals on the domain at X3 due to f\ and fi, respectively. 
Then T\\ and 7) 3 at the same crack front location may be expressed as 


4(E)1 = Ti * s 'i5 

4(* 3 )j Ts Tv 


7 a) (x 3 ) Ts / , 

— — ) 

f\ S 33 

7 (2) (E) 

/a 


The extension strain can be determined independently from finite element results. 



34 


For a monoclinic material about X3=0, s' 5 = 0 . Since f\ and fi, are arbitrary, they 


may be chosen as /j = f 3 = 1 without loss of generality. Therefore for this kind of 

material, at any given crack front location X 3 , T\\ and T\$ can be determined solely by / ll 
and P\ respectively. Eq.(4.7) then can be de-coupled as 


T u (x 3 ) = — 

y, 


r l> (x 3 ) — -r..(x.) 
s 33 


(4.8a) 


ri 3 (x 3 ) 



(4.8b) 


In modeling three-dimensional cracks along a given location X3, the I 33 
component is also induced by the extension strain along the crack front. Thus I 33 

can be evaluated as [37] 

TJx 3 ) =-[e }i (x 3 )-(s n T u + s }S T n )]. (4.8c) 

s 33 


For isotropic materials, 

v .. _ 1 „ _ a. „/ _ 1-1/2 „/ 2 (l + v) 

g ’ S 33 ~ g ’ S 35 ~ ® ’ N 1 — g ’ S 55 g 


Then Eq.(4.8) reduces to 

T i 1 0 3 ) = — [/ (1) (x 3 ) + ve n (x 3 )]; 

1 — V 


TJX,)= 2t^v) im{X ' y ' 

T 33 (x 3 ) — Ee 33 {x 3 ) + vTj^Xj) . 


(4.9) 


(4.10a) 

(4.10b) 

(4.10c) 


Note that Eq.(4. 10a) and Eq.(4. 10c) have been derived by Nakamura and Parks [28]. 
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5 Models 
5.1 Plates 

Due to loading and geometry symmetry, one-eighth of a through-thickness center- 
cracked plate is modeled by finite elements. The entire plate, shown in Figure 5.1, has a 
total length of 2/ = 80 in. , a total width of 2H’ = 40in. , and a total thickness of 
It = 0.33 in. The origin of the global Cartesian coordinate system is located at the center 
of the plate. The X-axis is parallel to the crack flank surfaces and the 7-axis is orthogonal 
to X and to the crack flank. The Z-axis is normal to the X-Y plane. A uniform 
displacement iu equivalent to a strain value of 0.1% is prescribed on the far ends at 
Y = ±1 . For the finite element model, symmetric boundary conditions are imposed on the 
planes of X = 0 , 7 = 0 (a < X < w ), and Z = 0 . The displacement loading is also applied 
on the surface at 7 = / . A local coordinate system is defined on the crack front at the 
centerline of thickness. The xi-axis is perpendicular to the crack front and coincident with 
the global X-axis. The x' 2 -axis is also normal to the crack front but parallel to the global 7- 
axis, and the xvaxis lies on the crack front and is parallel to global Z-axis as well, as 
shown in Figure 5.2. 

The finite element model is generated and solved by AN SYS, a general-purpose 
finite element code. A typical mesh of the finite element model is shown in Figure 5.3. It 
includes large brick-shaped elements far away from the crack front and small fan-shaped 
elements near the crack front. The elements attached on the crack front are 15 -node 
wedge-shaped elements. The radial size of these elements, eo, is always equal to or less 
than 1% of the characteristic length l c , defined as l c = min[<7, w - a, t] , i.e. e 0 < 0.01 xl c , 

as shown in Figure 5.4. Other than these wedge-shaped elements, the crack front region 
consists of 12 rings, each of which contains 12 elements. The width of each ring 
(equivalent to the radial element size) increases as the ring moves farther away from the 
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crack tip. The largest radial element size at Ring #12 is 26 times the smallest one at Ring 
#1, and its outer radius is kept at 100 times of eo, as shown in Figure 5.5. The half 
thickness is divided into 20 element layers, in which are 5 small layers near the free 
surface, 8 large layers near the centerline, and 7 mid-size layers in between. The three 
element sizes of each of these layers, in terms of the half thickness t, are 0.037, 0.057, and 
0.06257, respectively. Figure 5.6 shows the mesh refinement and element sizes near the 
crack front and through the thickness. The overall mesh contains 5640 elements and 
25183 nodes with 75549 degrees of freedom. This mesh result showed appropriate 
convergence on the calculation of the fracture parameters. 

5.2 Stiffened Panels 

A stiffened panel based on the dimensions of the plate mentioned in Section 5.1 
has five stiffeners attached longitudinally on one side of the plate. This structure is 
designed by Boeing for the all-composite wing skin in a commercial aircraft [29]. A 
center crack of length 2 a cuts through the central stiffener and the panel. Figure 5.7 
shows the configuration of the panel as well as the detail dimensions of the stiffener. Due 
to the presence of stiffeners and geometry symmetry, one-fourth of the entire panel is 
modeled by finite elements, as shown in Figure 5.8. Both of the global and local 
Cartesian coordinate systems are defined in the same way they are defined in a plate. 
Flence the symmetric boundary conditions are imposed on the planes of X = 0 and 7 = 0 
( a<X<w ). To prevent a free-body motion in Z-direction, a constraint in the Z- 
direction is also imposed on the node at (0, 0, -7), i.e. the point at the center of the panel's 
back side. Here 7 is defined as one half of the panel thickness that is not including the 
stiffener portion. A uniform displacement equivalent to a strain value of 0.1% is 
prescribed on the far end at Y = 1 . 

Since the crack is not expected to propagate across the “second stiffeners”, which 
are the stiffeners next to the central one, the crack front will only extend to near the edge 
of the second stiffeners. Under this assumption, the definition of the crack aspect ratio 
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will be different from what is defined in unstiffened plates. Instead of a/w, the crack 
aspect ratio in stiffened panels is defined as a/w’. Here a 'is the crack length calculated 
from the edge of the central stiffener to the crack front, and u’ ’ is the distance between 
edges of two adjacent stiffeners. Figure 5.9 shows an enlarged region along the crack 
surface with the designation of a, a \ and w’. 

The mesh pattern on the cracked panel portion of the finite element model is 
similar to that of the entire plate described in Section 5.1. The characteristic length l c , 
however, is defined as l c = min [<:/', w' - a',2t] . The radial size, eo, of those 15-node 
wedge-shaped elements attached on the crack front is set to be equal to 0.4% of the 
characteristic length, i.e. e 0 = 0.004/. . On the part of the finite element model containing 

uncracked panel and 3 stiffeners, the mesh is controlled to be as coarse as possible in 
order to reduce the overall number of elements. The full thickness of the panel, 2/, is 
divided into 30 element layers, in which are two symmetric parts that includes 15 layers 
each. For each half thickness, 5 large layers are near the centerline, 5 small layers are 
near the free surface, and 5 mid-size layers are in between. The three element sizes of 
each category of these layers, in terms of the half thickness / and from large to small, are 
0.04/, 0.07/, and 0.09/, respectively. Figure 5.10 shows the mesh refinement and element 
sizes near the crack front and through the thickness. The overall mesh contains 10539 
elements and 46780 nodes with 140340 degrees of freedom. 
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Figure 5.1 A through-thickness center-cracked plate subjected to a uniform far field 
displacement. 



Figure 5.2 One-eighth of the plate to be generated as a finite element model. 
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Figure 5. 6 Sizes of element layers in terms of the half thickness t. 
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Figure 5.9 Enlarged finite element mesh showing definition of the crack aspect ratio 
a/w’ (a/w’=0.1). 



Figure 5.10 Mesh refinements near the crack front of the stiffened panel (a Zw —0.1). 
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6 Results 


For the purpose of comparison, the stress intensity factors and 7-stresses are 
normalized as dimensionless numbers. The normalized stress intensity factor has the form 

of K t = K J (o' r while the two normalized 7-stresses are defined as T u = T u /a , 

and 7 13 =T u /(j x , respectively. Here o „ is the average stress, computed from the total 
nodal force, F due to the far-field displacement u x , on the plane of Y = 1 (see Figure 5.1) 
divided by the cross-sectional area of the plate. For a cracked plate, = F/(wt), where 
w is the half width and / is the half thickness. The thickness of the plate is also 
normalized because the variations of the stress intensity factor and the 7-stresses over the 
thickness are to be investigated. The normalized thickness of an element is defined as 
/ = x 3 // , which indicates the normalized X 3 coordinate at the centroid of the element. 

To determine the distribution of the stress intensity factor and the 7-stresses over 
the thickness of a particular model, these parameters are first calculated from each of the 
three adjacent rings of elements very close to the crack front (Rings #2, #3, and #4 as 
shown in Figure 5.4) by the equivalent domain integral and interaction integral. Then an 
average value is earned out over three domains. 

6.1 Isotropic Plates 

A family of isotropic plates with different crack lengths and different thicknesses 
is examined. The material properties of Y oung’s modulus 7 = 10 Msi and Poisson’s ratio 
v = 0.3 are used in the modeling. To compare the fracture parameters for different 
variables, the stress intensity factor K\ and T\\ stress are retrieved from the element layer 
attached to the centerline of the thickness. As shown in Figure 5.6, this layer has a 
thickness of 0.0625/, and the X 3 coordinate of the centroid is located at 0.03125/. Hence 
and 7 n of different models will be retrieved and plotted at normalized thickness 
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t= 0.03 125. The 7 b stress is retrieved from the element layer approximately at the 
quarter thickness of a plate, or at the thickness between 0.50t and 0.55 1, i.e.. Layer #9 as 
shown in Figure 5.6, and therefore t = 0.525 for this element. 

6.1.1 Crack Length 

First we will consider a plate of thickness 2t = 0.33in. (t/w = 0.00825 ) with 
various crack aspect ratios a/w = 0A~ 0.9 . To compare the effects of the thickness on 
fracture parameters for the same crack aspect ratios, the results from another set of plates 
with thickness of 2 1 = 10.24 in. ( // w = 0.256 ) are presented as well. 

Figure 6.1 shows the distribution of the normalized stress intensity factors K 1 
over half of the thickness for different crack aspect ratios. The distribution of the stress 
intensity factor is symmetric with respect to x?=0. K l is fairly constant over half of the 
thickness, except in the region near the free surface (approximately from 7 > 0.8 for 
a/w = 0A down to t > 0.6 for a/w = 0.9), where the stress intensity factor slightly 
decreases. Note that there exists a weaker corner singularity near the free surface. The 
normalized stress intensity factors at the center of thickness are retrieved and plotted 
against the a/w ratios in Figure 6.2. The results for t = 0.1 65 in. show that K l rises 
dramatically from 1.06 at the small crack of a/w = 0.1 to 2.72 as the alw ratio reaches 
0.9. Results from the other set of plates with thickness of t = 5. 12 in. (t/ w = 0.256 ) are 
also shown in the same figure, in which the normalized stress intensity factor increases 
from 1.00 to 2.64 as the a/w ratios goes from 0.1 to 0.9. It is slightly smaller than that of 
the thinner plate. 

The 2-D plane strain analytical solution of K l , also shown in Figure 6.2, has the 
following form in terms of the crack aspect ratio a/w [38]: 

Z : (— )= secf-(— )1 1 — 0.025C — f +0.06(— ) 4 . 
w \ 2 w JL w w 


( 6 . 1 ) 
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When the normalized stress intensity factor obtained from 3-D analysis is compared with 
the 2-D solution, i.e. Eq.(6.1), deviation within 4% on the average from the 2-D solution 
is found for the thicker plate; more than 5% for the thinner plate. 

Figure 6.3 shows the distribution of the normalized T\\ stresses over half of the 
thickness in the thinner plate (t = 0.165 in. ) for different crack aspect ratios. The T\\ 
stress distribution is symmetric with respect to X3=0. For all a/w ratios, the variation of 
T u over half of the thickness is relatively small except in the region near the free surface 
(approximately 7 >0.85), where the corner singularity would cause the T\\ stress to 
diverge. For most curves (a/w <0.7 ), the near-constant part of T n stays near -1.0, 
which is the theoretical 2-D solution. In Figure 6.4 the curve of the thinner plate 
(t/w = 0.00825 ) confirms this trend as the magnitude of the normalized T\\ stresses stays 
between -0.8 and -1.0 before the crack aspect ratio reaches 0.7. Beyond that the 
normalized Tu stress increases (in magnitude) rapidly to -3.77 as the a/w ratio reaches 
0.9. The other curve representing the thicker plate shows a similar trend, but the 
magnitude of T u at the center increases from -1.02 at the smallest crack, a/w = 0.1 , to 
- 4.87 at the largest crack, a/w = 0.9 . 

The 2-D solution of T u , also shown in Figure 6.4, can be obtained from the stress 
biaxiality ratio B and normalized stress intensity factors K t : 

T n =BK t , (6.2) 

where B is a function of the crack aspect ratio a/w [39]: 

£(— ) = - l + 0.085(— ) . (6.3) 

w i w 

Results from the thicker plate are closer to the 2-D solution than those of the thinner 
plate, although the closeness is only maintained until a/w is near 0.7. This indicates that, 
at a small crack, 7] , of the thicker plate is closer to the plane strain condition than that of 
the thinner plate. 
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In Figure 6.5, the distribution of the normalized Tn stresses over half of the 
thickness of the thinner plate for different a/w ratios does not have a constant portion 
shown in the distribution of K T and T n . Note that the distribution is anti-symmetric with 
respect to X3=0. Instead, the magnitude of the normalized Tn stresses rises monotonically 
from zero at the center of thickness to a much larger value near the free surface, and 
ultimately back to zero to satisfy the stress-free condition at the free surface. The larger 
the crack length is, the more sharply the 7j 3 stress results. As the normalized Tn stresses 

at t = 0.5 are retrieved for all a/w ratios, as shown in Figure 6.6, it is observed for the 
thinner plate that the magnitude of the normalized Tn stresses increases from 
approximately -0.35 for a/w = 0.\ to -2.6 for a/w = 0.9. For the case of the thicker 
plate, the trend persists. The overall variation of T u , however, is much smaller than that 
of the thinner plate. The absolute value of the normalized Tn stresses, which goes from 
approximately -0.07 for a/w = 0.1 up to near -0.6 for a/w = 0.9, is relatively smaller. 
Note that T n = 0 for the 2-D solution. 

6.1.2 Plate Thickness 

Another set of plates with different thicknesses is analyzed. The crack aspect ratio 
of these plates is kept at a/ w = 0.1. The finite element mesh of each model is also kept 
the same; i.e., the overall number of elements is unchanged for each finite element model. 
Only the dimensions of the elements are slightly changed due to thickness change. The 
thickness ranges from a thinner plate of 2/ = 0.33 m. (t/w = 0.00825 ) to a very thick 
plate of 2 1 = 20.48 in. (t/w = 0.512 ). 

Figure 6.7 shows the variation of normalized stress intensity factors over half of 
the thickness for plates with different t/w ratios. It is observed that, for relatively thin 
plates ( t/w < 0.064), K T tends to decrease as I goes away from the center of the plate. In 
contrast, for relatively thick plates, K 1 tends to increase in the thickness beyond the 
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region of t >0.3. The variation of K 1 for all t/w ratios at the center of thickness is 
plotted in Figure 6.8. For very thin plates (approximately t/w < 0.03), the normalized K\ 
almost stays constant near 1.06. As the thickness increases up to t/w > 0.5 , K l gradually 
decreases to near 1.0, which is the 2-D solution for very small crack length, i.e., 
a/ w — > 0 . 

The normalized T\\ stresses over half of the thickness for plates with different t/w 
ratios are plotted in Figure 6.9. The normalized T\\ stresses are fairly constant in the 
region of t < 0.75 for all cases. Beyond that, T n in a thinner plate deviates much more 
than it would in a relatively thick plate. These stresses, or the normalized T\\ stresses at 
center of the plate thickness, are extracted and plotted against the t/w ratios in logarithmic 
scale in Figure 6.10. The stresses, as shown also in Figure 6.9, increase gradually in 
magnitude from approximately -0.79 in a very thin plate ( t/w = 0.004 ) to -1 .01 in a very 
thick plate (t/w = 0.512). 

In Figure 6.11, the distribution of the normalized Tn stresses over half of the 
thickness rises in magnitude from zero at 1 = 0 to a larger value near the free surface. 
This trend is the same as that in Figure 6.5. In general the normalized Tn stresses 
decrease in terms of the magnitude as the plate thickness increases. The trend is observed 
for the normalized Tn stresses at t = 0.5 of each plate in Figure 6.12, in which the 
magnitude of T u decreases from approximately -0.5 for a very thin plate (t/w = 0.004 ) 
to near -0.02 for a very thick plate ( t/w = 0.5 1 2 ). 

6.2 Orthotropic Plates 

The plates of same dimensions are analyzed again with the properties of a 
composite material. The composite is AS4 carbon warp-knit fabric design by Boeing for 
the all-composite wing skin in a commercial aircraft [29], It has the overall orthotropic 
material properties that are listed in Table 6.1. The orientation of the material coordinates 
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Table 6. 1 Material properties of the orthotropic plate. 


E x = 5.162 Msi 

G yz = 0.64 Msi 

= 0.22 

E y =11.773 Msi 

Gjz = 0.57 Mpsi 

v xz = 0-29 

E z =1.53 Msi 

Gw = 2.479 Msi 

v M = 0.1758 


is referred to the global coordinate system shown in Figure 5.1. 

6.2.1 Crack Length 

The distribution of the normalized stress intensity factors K l through half of the 
thickness for various crack aspect ratios is shown in Figure 6.13. Each curve of the five 
different a/w ratios is relatively flat, meaning that K l is almost constant over the 
thickness. In fact, the amount of fluctuation of K l is less than 3% and occurs only in the 
region where t > 0.85. In comparison to Figures 6.13 and 6.1, it is observed that K l of 
the orthotropic material is more constantly distributed than that of the isotropic material. 
The magnitude for each a/w ratio for the orthotropic material is also slightly (about 6%) 
smaller. The normalized stress intensity factors at the center of thickness increases from 
1.00 to 2.54 as the a/w ratio increases from 0.1 to 0.9, as shown in Figure 6.14. Results 
from the thicker plate (t/w = 0.256 ) are also shown in the same figure, in which K 1 
ranges from 1.00 to about 2.52. Figure 6.14 shows that the two curves are almost 
identical, indicating that K 1 is less sensitive to the change of plate thickness in the 
orthotropic material. 

Figure 6.15 shows the distribution of the normalized T\\ stresses through half of 
the thickness for various crack lengths in the thin plate (t = 0.1 65 in. ). Similar to that of 
the isotropic material in Figure 6.3, the variation of T n is not significant except near the 
free surface (approximately t > 0.85 ). For the small crack case (a/w = 0.1 ), the flat part 
of the T u curve stays around -0.660, which is nearly equal to the material anisotropy 
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ratio in the 2-D solution [29], - ^js' 22 /s' u = -0.668. In Figure 6.16, the magnitude of T n 
at the center of the thickness increases gradually from -0.66 to -3.22 as the a/w ratio 
increases from 0.1 to 0.9, for the case of the thinner plate. The curve of the thicker plate 
(t/w = 0.256 ) is almost identical to that of the thinner plate, as it ranges between -0.67 
and -3.24. For the thicker plates in both materials, T u in the orthotropic material is about 
two-third of T u in the isotropic material (see also Figure 6.4 ). The 2-D solutions of both 
material types are about in the same ratio, i.e., -0.668 and -1.0. 

The trend of T u distribution over half of the thickness for various crack lengths in 

Figure 6.17 is similar to that of the isotropic material in Figure 6.5, except the magnitude 
for the orthotropic material is much smaller than that of the isotropic cases. In Figure 
6.18, normalized Tn stresses of the thinner plate at quarter thickness (t~ 0.5) range from 
about -3.94 xl0“ 2 to - 2.87x10“' as the a/w ratio increases from 0.1 to 0.9. T n of the 

thicker plate ranges approximately between -7.31xl0" 3 and -5.86xl0“ 2 . Both curves 
retain the trend that appeared in Figure 6.6. The magnitude of T u for the orthofropic 
material, however, is one order less than that in the isotropic material. 

6.2.2 Plate Thickness 

As the crack aspect ratio a/w is kept at 0.1, various plate thicknesses are adopted 
in the finite element model for analysis. The overall mesh configuration and total number 
of elements are unchanged, although the element size will be slightly changed due to 
thickness change. The half thickness varies from a very thin plate of 0.04in. 
(t/w = 0.002 ) to a very thick plate of 10.24in. (t/w = 0.512 ). 

The distribution of the normalized stress intensity factors through half of the 
thickness for plates with various t/w ratios is shown in Figure 6.19. The overall variation 
of K l for each thickness is relatively small (< 3%). A similar trend to the isotropic cases 
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is that K { tends to decrease in the region of t > 0.75 for thin plates ( t/w < 0.064 ) and 
increase in the region of t > 0.6 for thick plates. In the region of t < 0.5 , K 1 of each 
case stays very close to 1.0, the 2-D solution for very small crack lengths. K 1 at the 
center of thickness for all t/w ratios is retrieved and plotted in Figure 6.20, in which all 
normalized stress intensity factors fall between 0.995 and 1.005. The overall difference is 
within 1%. Therefore, K Y is rather stable even if the plate thickness increases more than 
two orders of magnitude. 

Figure 6.21 shows the normalized T\\ stresses over half of the thickness for 
various t/w ratios. For each plate, T u is fairly constant in the region of t < 0.75 . In the 
region near free surface (t> 0.85 ), however, T u tends to diverge as also observed in the 
isotropic cases (see Figure 6.9). It is observed that approximately at t = 0.8, all five 
curves pass through a point where T u is near -0.669 that is close to the material 
anisotropy ratio, -0.668. T u at the center of the thickness for all t/w ratios is extracted 
and plotted in Figure 6.22. For the plates of t/w > 0.00825 , T n ranges between -0.658 
and -0.673, or within ±1.5% of the material anisotropy ratio. As the t/w ratio increases 
beyond 0.016, i.e., plates with moderate to thick thicknesses, T u falls within ±1% of the 
material anisotropy ratio. As the plate becomes thinner (t/w <0.006), however, the 
magnitude of T u decreases from near the material anisotropy ratio to about 90% of that 
at t/w = 0.002 . 

In Figure 6.23, a similar trend to the T u distribution of isotropic plates over half 
of the thickness for various t/w ratios is observed. The relative magnitude of T u for 
orthotropic plates is one order less than that in isotropic plates. Figure 6.24 shows Ij 3 at 
quarter thickness for all plates. As the plate thickness increases to the very thick case of 
t/w = 0.512 , the magnitude of T u decreases to near zero. All normalized 7% stresses in 
orthotropic plates are relatively small, compared to those in isotropic plates. 
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6.3 Stiffened Panels 

All stiffened panels have the orthotropic material properties as listed in Table 6.1 
and the dimensions as shown in Figure 5.7. Panels with different crack aspect ratios, 
a Zw ’, ranging from 0.1 to 0.9 are analyzed. To illustrate the effect of crack aspect ratios 
on normalized stress intensity factor and the ^-stresses. each parameter at the centerline 
of the thickness will be plotted against a/w’ ratios. The value of a parameter is calculated 
from the average values over element layers #15 and #16, which are attached on the 
centerline. The thickness of each layer is 0.09/. Therefore the centroid is at -0.045/ for an 
element in Layer #15 and 0.045/ for an element in Layer #16. 

The original stiffened panel has a thickness of 2/ = 0.33 in. (t/w = 0.00825 ). To 
compare the effects of panel thickness on the fracture parameters with the same crack 
length, the results from another set of thinner panels of t/w- 0.004 are presented as 
well. Note that the stiffener dimensions are fixed for all the normalized studies. 

The distribution of the normalized stress intensity factors K Y through the entire 
thickness of the cracked panels is shown in Figure 6.25. For each crack aspect ratio, the 
distribution of K { appears to be increasing almost linearly from the bottom (7 = -l) to 
the top ( / = 1 ) of the panel thickness. The trend clearly shows the bending effect for 
thinner stiffened panels. The slope of each K l curve decreases as the a’/w’ ratio 
increases. For instance, K 1 increases from about 0.8 to 1.9 for a'/w' = 0.1 but ranges 
only between 0.96 and 1.06 for a'/w' = 0.9. Figure 6.26 shows K 1 at center of the 
thickness for panels with various crack aspect ratios. The normalized stress intensity 
factor of the thicker stiffened panel (t/w = 0.00825 ) decreases from approximately 1.38 
to 1.05, as the a’/w’ ratio increases from 0.1 to 0.9. The thin stiffened panels 
(t/w = 0.004 ) generally have larger K 1 , which decreases from 1.67 to 1.05 for a’/w’ ratio 
from 0.1 to 0.9. The difference of two sets of panels, however, becomes smaller as the 
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crack length increases. Both panels almost have an equal K 1 of 1.05 at a'/w' = 0.9. It 
indicates that, as the crack length becomes large enough, the panel thickness is irrelevant 
to the normalized stress intensity factor. 

Figure 6.27 shows the distribution of the normalized T\\ stresses over the 
thickness of the panel for different a'/w' ratios. The magnitude of T n is generally in an 
increasing trend, although the distribution appears more linearly in larger crack lengths. 
Overall, the magnitude of T n decreases as the crack length increases. This can be 
observed in Figure 6.28 where the normalized Tn stresses are extracted and plotted 
against different a Zw’ ratios. Though all are negative, T u of the thicker panel increases 
from -0.91 at a'/w' = 0.1 to -0.52 at a'/w' = 0.9, while T u of the thinner panel ranges 
from -1.17 to -0.41. For smaller cracks ( a'/w' < 0.7 ), the normalized T\\ stress of the 
thinner panel is larger (in magnitude) than that of the thicker panel. For larger cracks 
(a'/w' >0.8), however, the trend reverses as the thinner panel has a smaller T n (in 
magnitude). 

In Figure 6.29, the distribution of T u through the thickness appears in a nearly 
anti-symmetric manner with respect to the thickness centerline, t = 0 . The shape of each 
T n curve is similar enough that it seems all curves are shifting within a range 
approximately equal to 0.5. The overall variation of T u for all c//w’ ratios ranges between 
-0.1 and -1.5. T u at center of the thickness for each crack length of both sets of panels is 
retrieved and plotted in Figure 6.30. Both curves form part of a parabola, respectively. 
For the thicker panels (t/w = 0.00825 ), T u decreases from -0.93 at a'/w' = 0.1 to -0.95 
at a'/w' = 0.3 near the lowest point of the curve, then increases to about -0.58 at 
a'/w' = 0.9. For the thinner panels, the normalized Tn stress swings from -1.86 at 
a'/w' = 0.1 to -1.22 at a'/w' = 0.9. In between the lowest point is near a'/w' = 0.4 
where T u is approximately -2.16. Figure 6.30 shows the thinner panel has larger T u (in 
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magnitude) for different crack lengths. And the maximum amplitude of T u for a panel of 
fixed thickness occurs at a moderate crack aspect ratio. 




Figure 6.2 Normalized stress intensity factors at center of the thickness (xft=0) for 
isotropic plates with Wo different thicknesses. 
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Figure 6.3 Normalized T u stresses through half of the thickness for isotropic plates of 


t=0.165 in. (t/w=0.00825) with various a/w ratios. 



Figure 6.4 Normalized T\\ stresses at center of the thickness (xft=0) for isotropic 
plates with Wo different thicknesses. 
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Figure 6.5 Normalized Tn stresses through half of the thickness for isotropic plates of 
t=0.165 in. (t/w=0. 00825) with various a/w ratios. 



Figure 6.6 Normalized stresses at quarter of the thickness (xft=0.5) for isotropic 
plates with Wo different thicknesses. 
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Figure 6 . 7 Normalized stress intensity factors through half of the thickness for 


isotropic plates of a/w=0. 1 with various t/w ratios. 



Figure 6.8 Normalized stress intensity factors at center of the thickness (xft=0) for 
isotropic plates of a/w=0. 1. 
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Figure 6.9 Normalized T\\ stresses through half of the thickness for isotropic plates of 
a/w=0. 1 with various t/w ratios. 



Figure 6.10 Normalized T\\ stresses at center of the thickness (xft=0) for isotropic 
plates of a/w=0.1. 









Figure 6.14 Normalized stress intensity factors at center of the thickness (xft=0) for 
orthotropic plates with Wo different thicknesses. 
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Figure 6.15 Normalized T\\ stresses through half of the thickness for orthotropic plates 


of t=0.165 in. (t/w=0.00825) with various a/w ratios. 



Figure 6.16 Normalized T\\ stresses at center of the thickness (xft=0) for orthotropic 
plates with Wo different thicknesses. 





62 



Figure 6.17 Normalized Tn stresses through half of the thickness for orthotropic plates 


of t=0.165 in. (t/w=0.00825) with various a/w ratios. 



Figure 6. 18 Normalized Tn stresses at quarter of the thickness (xft=0.5) for orthotropic 
plates with Wo different thicknesses. 




63 



Figure 6.19 Normalized stress intensity factors through half of the thickness for 


orthotropic plates of a/w=0. 1 with various t/w ratios. 



Figure 6.20 Normalized stress intensity factors at center of the thickness (xft=0) for 
orthotropic plates of a/w=0. 1. 
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Figure. 6.21 Normalized T\\ stresses through half of the thickness for orthotropic plates 
of a/w=0. 1 with various t/w ratios. 



Figure 6.22 Normalized T\\ stresses at center of the thickness (xft=0) for orthotropic 
plates of a/w=0.1. 
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Figure 6.25 Normalized stress intensity factors through the thickness for orthotropic 
stiffened panels oft=0.165 in. (t/w=0.00825) with various a/w’ ratios. 



Figure 6.26 Normalized stress intensity factors at center of the thickness (xft=0) for 
orthotropic stiffened panels with various a/w ’ratios for two thicknesses. 
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Figure 6.27 Normalized T\\ stresses through the thickness for orthotropic stiffened 
panels of t=0.165 in. (t/w=0. 0082 5 J with various a Zw’ ratios. 



Figure 6.28 Normalized T\\ stresses at center of the thickness (xft=0 ) for orthotropic 
stiffened panels with various a/w’ ratios for two thicknesses. 
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Figure 6.29 Normalized Tn stresses through the thickness for orthotropic stiffened 
panels of t=0.165 in. (t/w=0.00825) with various a Zw’ ratios. 



Figure 6.30 Normalized 7 b stresses at center of the thickness (xft=0 ) for orthotropic 
stiffened panels with various a Zw’ ratios for two thicknesses. 
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7 Summary and Conclusions 


The equivalent domain integral method is used to determine the point-wise stress 
intensity factor, or the distribution of it along a three-dimensional crack front. This 
method is modified from ./-integral with the use of a weighting function, i.e., the s- 
function. In the practice of numerical computation, the finite element shape function is 
applied as the ^-function. 

In addition to the stress intensity factor, a second fracture parameter, namely the 
T-stress, is determined in order to better characterize the fracture behavior. The 
components of T-stress can be obtained by the evaluation of the interaction integral. The 
interaction integral is derived from the equivalent domain integral method with the 
concept of an auxiliary field. The auxiliary field is the solution corresponding to uniform 
force acting on the crack front. The Stroh formalism is used to derive the auxiliary stress 
and displacement fields associated with different T-stress components, such as T\\ and 
T]3, in the monoclinic material with the plane of symmetry at X3=0. 

The finite element models are made on two sets of through-thickness center- 
cracked plates with isotropic and orthotropic material properties, respectively. Similar 
finite element meshes are generated on another set of orthotropic composite panels with 
stiffeners. All of these structures are under the mode-I uniform displacement loading. 

For plates, the stress intensity factor will increase as the crack length increases in 
plates with the same dimension. The effect of material properties on the stress intensity 
factors at the center of the plate is not significant and all of the stress intensity factors are 
close to the corresponding 2 -D solutions, especially for relatively thick plates. For plates 
of different thicknesses with the same crack length, the stress intensity factor for isotropic 
plates decreases as the thickness increases while for orthotropic plates the trend is 
relatively insensitive. For panels with stiffeners, however, the stress intensity factor 
decreases as the crack length in the panel increases. Similar to the unstiffened plates, a 
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thinner stiffened panel has a larger stress intensity factor for a given crack length. The 
exception occurs for the case of relatively larger crack lengths, e.g., a'/w' = 0.9, where 
two sets of thin panels of different thicknesses have nearly identical values of stress 
intensity factor. 

The T\\ stresses in all cases of center-cracked plates and panels are compressive. 
Based on the stress field near the crack tip prior to fracture initiation from a couple tests 
[29], the crack may not have tendency to turn. This is in contradiction to the experimental 
observation that the crack did turn near the stiffener. A possible reason may be due to 
other interacting failure modes, such as delamination and disbond with stitched interface, 
that caused the turning. Among plates, the magnitude of T\\ increases when the crack 
length increases, regardless of the material properties or plate thickness. The thickness 
effect in plates with a small length of crack is more profound for isotropic ones, in which 
the T\\ stress increases (in magnitude) as the thickness increases. Although in general the 
orthotropic plates have the same trend, the changes of T\\ are relatively small for thicker 
plates. The trends observed in plates are reversed in stiffened panels, where T\\ decreases 
(in magnitude) when a crack length in a panel increases. At small crack lengths, T u in a 
thin panel also has a larger magnitude then Tu in a thick panel. 

The Tu stress in plates of both materials increases (in magnitude) when the crack 
length increases, but a larger magnitude is found in isotropic plates. Tu, however, 
decreases in magnitude as the plate thickness increases. For very thick plates, where the 
plane strain condition dominates near the middle of plate thickness, Tu approaches zero. 
This is also true in two-dimensional analysis, in which Tu is the only effective 7-stress 
component. In stiffened panels, the magnitude of Tu is generally smaller for large cracks. 
The largest magnitude of Tu , however, is not found in panels with a smallest crack, but 
with a somewhat larger crack. 

This research may be extended in the future according to the following aspects. 

(1) The loading condition may be extended from uniaxial tension to others, such as pure 
bending or biaxial loading. The change of the loading conditions, however, may incur 



71 


mixed-mode effects. If the asymptotic stress field near the crack front includes mixed 
modes, the mode-II and mode-III effects must be taken into account for the equivalent 
domain integral formulation. The calculation of the stress intensity factors will also 
require the consideration of the auxiliary fields [21], 

(2) The effects of mixed mode may come from the existence of a kinked crack [40]. 
However, a fracture plate with the kinked crack requires an even larger finite element 
model to simulate because it usually is not considered as a symmetric geometry. 

(3) Since the stress intensity factors and the T-stresses are related to the material 
properties, a composite material may be designed such that a certain level of the 
desired parameters can be obtained. 
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A. ANSYS Program 


The following ANSYS program is to generate a finite element model of one- 
fourth of a stiffened panel with a center crack. The structural dimensions are described in 
Chapter 5, and the material properties are listed in Table 6.1. The program will output 
plain text files which includes data of displacements, strains, stresses, nodal coordinates, 
and other parameters to be used later in a FORTRAN program. This program is tested on 
ANSYS/ResearchFS, the product combination of ANSYS Revision 5.5.1 licensed to the 
North Carolina Supercomputing Center. This program can be modified to generate a 
finite element model for the simpler center-cracked plate. 

ANSYS commands are not case-sensitive, i.e., a command in uppercase is 
identical to a command in lowercase. For the purpose of clarity, however, all commands 
and associated arguments are listed with uppercase characters, except the user-defined 
parameters that use lowercase characters. 
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/BATCH 

/UNITS, BIN ! British system of units using inches 

! ================ BEGINNING of Parameters Definition ================= 

! Material Properties 

*DIM,mat , ARRAY, 3 , 3 , 1 ! array to store material properties 

mat (1,1) = 5. 162e6, 11. 773e6, 1 .530e6 ! Ex,Ey,Ez 

mat (1,2) = 0.22,0.29,5.162/11.773*0.401 ! NUyz , NUxz ,NUxy 

mat (1,3) = 0 . 640e6 ,0.570e6,2.479e6 ! Gyz,Gxz,Gxy 

! Parameters for Loading Conditions 

appstn = 0.001 ! applied far-field uniform strain 


! Geometric Parameters 
*DIM, cir , ARRAY, 5,1,1 


aow 

= 

0. 1 

plz 

= 

0.33 

plx 

= 

20 

ply 

= 

40 

tst 

= 

0.22 

hst 

= 

2 . 3 

wst 

= 

1.6+tst/2 

wsl 

= 

wst-tst 

dxc 

= 

8 

dxe 

= 

dxc-2*wst 

dsp 

= 

ply*appstn 

clp 

= 

aow* dxe 

ckl 

= 

wst+clp 

ctf 

= 

dxe -clp 


! array of dimensions of the fan- type area 
! a/w ratio 

! full length of the panel in Z-direction (thickness) 
! 1/2 length of the panel in X-direction (w) 

1 1/2 length of the panel in Y-direction (1) 

! 1/2 thickness of the stiffener 

! height of the stiffener 

1 1/2 width of the stiffener (aO) 

1 1/2 width of the base of the stiffener 
1 distance between the center of 2 stiffeners 
1 distance between the edge of 2 stiffeners 
1 applied far-field displacement 
1 crack length in panel portion (a' ) 

! 1/2 crack length (a=a0+a' ) 

1 panel length in front of crack tip 


! Paramet 
*DIM,neci 
*DIM, srrl 
*DIM , zkpn 
*DIM , zntk 
*DIM,nezi 
neci ( 1 ) = 
srrl ( 1 ) = 
zntk(l) = 
ne z i ( 1 ) = 
narc = 6 
angc = 90 
nezz = ne 


ers for Mesh 
, ARRAY ,3,1,1 
, ARRAY ,3,1,1 
, ARRAY ,4,1,1 
, ARRAY ,6,1,1 
, ARRAY ,6,1,1 
1 , 12,1 
1,26,1 

0.20,0.35,0. 
5, 5, 5, 5, 5, 5 

/narc 

zi (1) +nezi (2) 


nerr = neci (2 ) 
neyf = 6 


sryf 


32 


Control 

1 mesh-control parameters over the fan 
1 spacing ratios for radial lines 
! keypoint numbers at Z-direction interfaces 
! normalized thicknesses for different spacing 
! number of elements in each region 


45,0.45,0.35,0.20 

! number of arcs for a 90 -degree span 
! arc angle per element 

+nezi(3) ! number of element layers along 1/2 

thickness 

! number of elements in r-direction 
! number of solid model divisions in 
Y-direction of far-field portion 
! spacing ratio in Y-direction of far-field portion 


! Define spacing ratios in X-direction 

*IF, aow, LE, 0.3, THEN 
neci = 2 
srl = 3 
nec2 = 4 
sr2 =10 

*ELSEIF, aow , GE , 0 . 7 , THEN 
neci = 4 
srl = 10 
nec2 = 2 
sr2 = 3 

*ELSE ! 0.3 < a'/w' < 0.7 
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necl = 3 
srl = 9 
nec2 = 3 
sr2 = 9 
*ENDIF 


! Define critical length 
*IF, aow, LE, 0.5, THEN 
* IF , plz , LE , clp , THEN 
lcr = plz 
*ELSE 
lcr = clp 
*ENDIF 
*ELSE 

* IF, plz , LE, ctf , THEN 
lcr = plz 
*ELSE 
lcr = ctf 
*ENDIF 
*ENDIF 


a' /w' <= 0.5 
(2t<=a' ) 

critical length 
(2t>a' ) 

critical length 

a' /w' > 0.5 
(2t<=w' -a' ) 
critical length 
(2t>w' -a' ) 
critical length 


min [a' , w' -a' 

,2t 

min [a' , w' -a' 

, 2 1 


min [a' , w' -a' 

, 2 1 

min [a' , w' -a' 

, 2 1 


rac = 0.004*lcr ! radial size of the smallest element (eO) 

cir(l) = 0 , rac, 100 *rac , 150*rac 

dyf = ply-cir(4) ! length in Y-direction 

tw = (plz/2)/plx 1 t/w ratio 

i END of Parameters Definition 


/TITLE, 1/4 stiffened panel , a' /w' =%aow% , t/w=%tw% , lcr=%lcr% , e0=%rac% ; C. Lin 


/PREP7 

ET, 1, 95, , , , , 1 

i 

KEYOPT, 1,11,1 
ET,2 , 93 

MP , EX, 1 , mat (1,1) 

MP , E Y , 1 , mat (2,1) 

MP , EZ , 1 , mat (3,1) 

MP , PRXY , 1 , mat (3,2) 

MP, PRYZ, 1 ,mat (1,2) 

MP, PRXZ, 1 ,mat (2 ,2) 

MP , GXY , 1 , mat (3,3) 

MP , GYZ , 1 , mat (1,3) 

MP , GXZ , 1 , mat (2,3) 

! Solid Model Generation: 

! Generate fan- type areas over a span of 90 degrees 
LOCAL, 14 , 0 , ckl , 0 , -plz/2 ! local Cartisian c.s. #14 at crack tip 
WPCSYS,1,14 ! set working plane at local c.s. #14 

*D0,i,l,narc ! generate the innermost ring. Ring #0 

PCIRC, cir(l) , cir (2) , (i-1) *angc , i*angc 
*ENDDO 

PCIRC, cir (2) , cir (3) , 0, 90 ! annular area to be mapped meshed later 

RECTNG, 0, cir (4) , 0, cir (4) 

ARSYM, X, ALL 

BTOL,rac/100 ! set Boolean operation tolerance 

AOVLAP,ALL ! create areas on the x-y plane 

BTOL , DEFA 


BEGINNING of Preprocessing Phase ================== 

! element type: SOLID95, the 20 -node 3-D structural 
element with solution output at integration points 
! 2x2x2 reduced integration option 
! SHELL93 , 8-node plate element 
1 assign material properties 


! Set line divisions on all circumferencial and radial lines 
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LOCAL, 15 , 1 , ckl , 0 , 0 ! local cylindrical c.s. #15 at crack tip 

LSEL , S , LOC , X , ci r (2 ) 

LESIZE,ALL, , , 1 ! set divisions on circumf erencial lines 

*DO, i, 1,3 

LSEL, S, LOC, X, (cir (i) +cir (i+1) ) /2 

LESIZE, ALL, , , neci (i) , srrl ( i) ! set divisions on radial lines 

*ENDDO 

! Concatenate lines around the fan-type area for mapped mesh later 
*DO , i, 1,2 

LSEL, S, LOC, X,cir(2) 

LSEL, R, LOC, Y, (i-1) *90, i*90 
LCCAT , ALL 

LSEL, S, LOC, X, cir (4) /2*SQRT (5) 

LSEL, R, LOC, Y, ( i - 1 ) * 90 , i* 90 
LESIZE, ALL, , ,narc/2 , 1 
LCCAT , ALL 
*ENDDO 
LSEL, ALL 

NUMCMP,ALL ! compress all numbering 

! Mesh areas on X-Y plane 

/NERR,3 ! turn off warning messages 

ASEL , S , LOC , X, 0 , cir (2) 

AATT ,1,1,2 

MSHAPE,1,2D ! mesh areas containing crack tip with triangle elements 

AMESH,ALL ! areas meshed with SHELL93 elements 

ASEL, INVE 
AATT ,1,1,2 

MSHAPE,0,2D ! mesh other areas with quadrilateral elements 

AMESH,ALL ! areas meshed with SHELL93 elements 

LSEL , S , SPACE ,, 0 ! select concatenated lines to be deleted 

LDELE,ALL ! delete concatenated lines 

! Generate the finite element model along full panel thickness 
ASEL, ALL 

TYPE,1 ! designate element type as SOLID95 

CSYS,0 ! set active c.s. as global Cartesian 

zi = -plz/2 

*D0,i,l,6 ! generate SOLID95 elements by extending selected areas 

ESIZE , , nezi (i) 

VEXT, ALL, , , , , zntk (i) *plz/2 

zi = zi+zntk (i) *plz/2 ! Z-coordinate of the areas to be offset 
ASEL, S, LOC, Z, zi 
*ENDD0 

ASEL, S, TYPE, ,2, , ,1 

ACLEAR , ALL ! delete all SHELL 9 3 elements 

! Construct the other portions of the panel 

CSYS,14 ! switch to the local Cartesian c.s. on the crack tip 

NUMSTR,KP, 1001 

NUMSTR, LINE, 2001 

KSEL, S,LOC, Z, 0 

KSEL , R , LOC , Y , 0 

KSEL, R, LOC, X, -cir(4) 

*GET, nkpl , KP, 0 , NUM,MAX ! nkpl : keypoint # of this keypoint 

KGEN, 2 , ALL, , , - (clp-cir (4 ) ) ! keypoint #1001 

L, nkpl , 1001 , necl , srl ! line #2001 

KSEL, S, LOC, Z, 0 

KSEL, R, LOC, Y, 0 
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KSEL , R , LOC , X , cir (4) 

*GET , nkp2 , KP , 0 , NUM , MAX ! nkp2 : keypoint # of this keypoint 

KGEN, 2 , ALL, , , dxe-clp-cir (4) ! keypoint #1002 

L, nkp2 , 1002 , nec2 , sr2 ! line #2002 

NUMSTR , DEFA 

ASEL , S , LOC, X, -cir (4) 

VDRAG, ALL, ,,,,, 2001 ! generate cracked portion of the panel 

ASEL, S, LOC, X,cir(4) 

VDRAG, ALL, ,,,,, 2002 ! generate uncracked portion of the panel 

! Generate the central stiffener 

CSYS , 0 ! switch back to global Cartesian c.s. 

ASEL, S, LOC, X,wst 
ESIZE, , 1 
VEXT , ALL , , , -wsl 
ASEL, S, LOC, X,tst 
ESIZE, , 1 
VEXT, ALL, , , -tst 
ASEL, S, LOC, Z,plz/2 
ASEL, R, LOC, X, 0, wst 
ESIZE, , 1 

VEXT, ALL, , , , , tst 
ASEL , S , LOC ,Z,plz/2+tst 
ASEL, R, LOC, X, 0, tst 
ESIZE, , 1 

VEXT, ALL, , , , , hst-tst 
! Generate the second stiffener 

VSEL , S , LOC, X, tst , wst , , 1 ! select volumes on the outer central 

! stiffener 

LOCAL, 11 , 0 , dxc/2 , 0 , 0 ! set local Cartesian c.s. #11 

VSYMM,X,ALL ! generate volumes on the inner second stiffener 

LOCAL, 12 , 0 , dxc, 0 , 0 ! set Cartesian c.s. #12 at the center of the 

! second stiffener 

ASEL, S, LOC, X, -tst 
ESIZE, , 1 

VEXT, ALL, ,, 2*tst ! create volumes in the base of the second stiffener 

ASEL , S , LOC ,Z,plz/2+tst 
ASEL, R, LOC, X, -tst, tst 
ESIZE, , 1 

VEXT , ALL ,,,,, hst -tst ! generate volumes in the height of the second 

! stiffener 

VSEL, S, LOC, X, -wst, -tst, , 1 

VGEN, 2 , ALL, , , wst+tst ! generate volumes in the outer part of the 

! second stiffener 

! Generate the outmost stiffener 

VSEL , S , LOC, X, -wst , wst ,, 1 ! select all volumes in the second stiffener 

LOCAL, 13 , 0 , 1 . 5*dxc , 0 , 0 ! set Cartesian c.s. #13 

VSYMM,X,ALL ! generate the entire outmost stiffener 

! Generate the uncracked panel and the outmost panel 
LOCAL, 16 , 0 , 2*dxc , 0 , 0 ! set Cartesian c.s. #16 

ASEL, S, LOC, X, -wst 
ASEL, R, LOC, Z, -plz/2,plz/2 

VEXT, ALL, , , -dxe ! generate the uncracked part of the panel 

ASEL, S, LOC, X, wst 
ASEL, R, LOC, Z, -plz/2,plz/2 
VEXT, ALL, , , dxe/ 2 


generate the outmost part of the panel 



ALLSEL , ALL , ALL 
NUMMRG , ALL , le-6 
NUMCMP , ALL 
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! compress all numbering 


! Construct the entire FE model by extending areas along Y-axis 


CSYS , 0 

NUMSTR,KP, 1001 
NUMSTR, LINE, 2001 
KSEL , S , LOC, X, 0 
KSEL , R , LOC , Z , 0 
KSEL, R, LOC, Y, cir (4) 
*GET , nkp3 , KP , 0 , NUM , MAX 
KGEN , 2 , ALL , , , , dyf 
L, nkp3 , 1001 , neyf , sryf 
ASEL , S , LOC, Y, cir (4 ) 
VDRAG, ALL, , , , , ,2001 


! nkpl : keypoint # of this keypoint 
! keypoint #1001 
! line #2001 

1 the entire finite element model completed 


ALLSEL, ALL, ALL 
NUMMRG, ALL, le-6 

NUMCMP, ALL ! compress all numbering 

/VIEW, 1 , 1 , 1 , 1 1 set viewing point for plots 

/ANGLE , 1 , - 9 0 , XM 
/AUTO , 1 
/TYPE, ALL, 2 


ASEL, S, LOC, Y, ply 
ASUM 

*GET , ay t , AREA , 0 , AREA 

/OUTPUT, par, dat 

* STATUS, ay t 

/OUTPUT 

ASEL, ALL 

FINISH 


1 calculate the total area on the far-field end 

! output data in <par.dat> file 
! "ayt": total area on the far-field end 


! finish generation of the finite element model 
END of Preprocessing Phase 


/SOLU 

ANTYPE, STATIC 

NSEL , S , LOC , X, 0 

D , ALL , UX , 0 

NSEL, R, LOC, Y, 0 

NSEL, R, LOC, Z, -plz/2 

D , ALL , UZ , 0 

NSEL, S, LOC, Y, ply 

D , ALL , UY , dsp 

NSEL, S, LOC, Y, 0 

NSEL , R , LOC , X , ckl , plx 

D , ALL , UY , 0 

NSEL, ALL 

EQSLV, PCG, le-8 

ERESX, NO 

SAVE 

SOLVE 

FINISH 


BEGINNING of Solution Phase ==================== 

1 static analysis 

1 set X-symmetry boundary conditions 

! set Z- constraints at the center of plate 
! apply the prescribed displacement at far end 

! set Y- symmetry boundary conditions 
! use the PCG solver 

! ouput integration point results to the nodes 
END of Solution Phase 


! ================= BEGINNING of Postprocessing Phase 

/POST1 

nezt = nezz*2 
*DIM , cmr ing , CHAR , 3 ,1,1 


! number of element layers along thickness 
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*DIM, zf , ARRAY, nezt+1 , 1 , 1 ! array to store Z-coordinate of every layer 

cmr ing ( 1 ) = ' ring2 ' , ' r ing3 ' , ' ring4 ' 

srr = srrl (2) ** ( 1/ (neci (2 ) - 1) ) 1 spacing ratio between adjacent 

! elements in r-direction 

esr = (cir (3) -cir (2) ) * (1-srr) / ( l-srr**neci (2) ) ! element size in r-direction 

zf(l) = -plz/2 ! Z-coordinate at bottom of the first layer 

kring = 3 ! number of rings to be recorded 

k = 1 

*D0,i,l,6 ! i: index for mesh-size region (= dimension of "zntk") 

eszi = zntk ( i) * (plz/2 ) /nezi ( i) ! element size (z) in each region 

*D0 , j , 1 , nezi (i) ! j: index for the number of element layer 

k = k+1 

zf (k) = zf(k-l)+eszi 
*ENDDO 
*ENDDO 


SET, 1 

NSEL , S , LOC, Y, ply 
FSUM 

*GET , f yt , FSUM , 0 , ITEM , FY 

/OUTPUT , par , dat , , APPEND 

* STATUS, fyt 

*STATUS , nezt 

*STATUS , plx 

♦STATUS, plz 

♦STATUS, mat, 1,3, 1,3 

♦STATUS, appstn 

♦STATUS, ckl 

♦STATUS, narc 

♦STATUS, kring 

♦STATUS, nerr 

♦STATUS, aow 

/OUTPUT 

NSEL, ALL 

CSYS, 15 

DSYS, 15 


! select nodes on the far end 
! calculate total nodal force on the far end 


I "fyt": total nodal force on the far end 
1 number of elements layers along thickness 


1 activate local cylindrical c.s. at crack tip 
! set display c.s. to #15 


! Write stresses 
! every elements 
/OUTPUT, sts, dat 
/OUTPUT 

/OUTPUT, stn, dat 
/OUTPUT 

/OUTPUT, dis, dat 
/OUTPUT 

/OUTPUT, node, dat 
/OUTPUT 


strains, displacements, and nodal coordinates of 
into files 

! <sts.dat> file to store stresses 
! <stn.dat> file to store strains 
! <dis.dat> file to store displacements 
! <node.dat> file to store node coordinates 


♦DO, k, 2 , kring+1 ! select from Ring #2 to Ring #(kring+l) 

NSEL, S, LOC, X, rac+esr* (1-srr** (k-1) ) / (1-srr) , rac+esr* (1-srr **k) / (1-srr) 
CM, cmring (k- 1) , NODE ! group all nodes in the domain to be integrated 

*DO,i,l,nezt ! i: index for layer in Z -direction 

NSEL, R, LOC, Z, zf (i) , zf (i+1) 

CM , nlayer , NODE 

♦DO, j , 1 , narc*2 ! j: index for element in the same layer 

NSEL, R, LOC, Y,angc* ( j - 1 ) ,angc*j 
ESLN,S,1 ! select an element 

/FORMAT, , E 

/HEADER , OFF , OFF , OFF , OFF , ON , ON 
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/OUTPUT, sts, dat , , APPEND 
*MSG, INFO, k, i , j 

Ring # %I , Segment # %I , Element # %I : 

PRESOL, S 
/OUTPUT 

/OUTPUT, stn, dat, , APPEND 
*MSG, INFO, k, i, j 

Ring # %I , Segment # %I , Element # %I : 

PRESOL , EPEL 
/OUTPUT 

/HEADER, OFF , OFF , OFF , OFF , OFF , OFF 
/OUTPUT, dis, dat, , APPEND 
*MSG, INFO, k, i , j 

Ring # %I , Segment # %I , Element # %I : 

PRNSOL , DOF 
/OUTPUT 

/OUTPUT , node , dat , , APPEND 
*MSG, INFO, k, i , j 

Ring # %I , Segment # %I , Element # %I : 

ELIST, ALL, ,,0,0 
NLIST, ALL, , , , Z, Y,X 
/OUTPUT 

CMSEL , S , nlayer 
*ENDDO 

CMSEL, S,cmring(k-1) 

*ENDDO 

*ENDDO 

! Write Strain_33 for every crack-front element. 

/OUTPUT, stn33 , dat ! <stn33.dat> file to store strains e_33 

/OUTPUT 

NSEL , S , LOC, X, 0 , rac 

CM, fan, NODE ! group all nodes in the fan- type domain 

/FORMAT, , E 

/HEADER , OFF , OFF , OFF , OFF , ON , ON 

*DO,i,l,nezt ! i: index for layer in Z-direction 

NSEL, R, LOC, Z, zf (i) , zf (i+1) 

CM, nlayer, NODE 

NSEL, R, LOC , X, 0 ! select 3 nodes on the crack front 

CM , nf ront , NODE 
CMSEL, S, nlayer 

*DO , j , 1 , narc*2 ! j: index for element in the same layer 

NSEL, R, LOC, Y,angc* (j-1) ,angc*j 
*IF, j , GT, 1, THEN 
CMSEL, A, nf ront 
*ENDIF 

ESLN,S,1 ! select an element 

/OUTPUT, stn33, dat, , APPEND 
*MSG, INFO, 0, i, j 

Ring # %I , Segment # %I , Element # %I : 

PRESOL, EPEL 
/OUTPUT 

CMSEL, S, nlayer 
* ENDDO 
CMSEL, S, fan 
* ENDDO 


/HEADER, DEFA 
NSEL, ALL 
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ESEL , ALL 
DSYS , 0 
CSYS , 0 
FINISH 

i END of Postprocessing Phase 

/EXIT 
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B. FORTRAN Program 


The following FORTRAN program is used to calculate the stress intensity factor 
and the T-stresses along a three-dimensional crack front in the stiffened panel. The crack 
is under mode-I loading (either a uniform displacement or a uniform stress) on the far 
ends of the structure. The material type can be isotropic or orthotropic. The general input 
of the program is from data files containing stresses, strains, displacements, nodal 
coordinates and other associated parameters that are generated from a finite element 
analysis run by ANSYS. 

The program is created and tested in the Microsoft Fortran PowerStation 4.0, a 
commercially available FORTRAN development software for Windows NT and 
Windows 95 operating systems. Although the Fortran PowerStation 4.0 incorporates all 
the features of Fortran 90, this program is intentionally written by the FORTRAN 77 
syntax. This arrangement would maintain the flexibility that the program can be also run 
on the platforms that still carry FORTRAN 77. However, the program employs some of 
the IMSL libraries, the integrated mathematical and statistical functions, which may not 
be available on every platform with the FORTRAN capability. 
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PROGRAM T_STRESS 

C This program calculates the stress intensity factor K_I and the T- stress 

C components of Til, T13 and T33 of a 3-D crack problem. 

C The capabilities of the program in calculation of the finite element 

C model are 50 element layers along the crack front, 24 elements per 180 
C degrees. 

C Definition of Selected Variables: 

C ADISP(i,j) Array for the auxiliary displacement field (for Til) within an 
C element. i=1..20. j=1..3. 

C ADISP13(i,j) Array for the auxiliary displacement field (for T13) within 

C an element. i=1..20. j=1..3. 

C APPSTN Far-field strain. 

C APPSTS Applied far-field stress. 

C ARC Arc length (in degree) of a single element 

C ASTRN(i,j) Array for the auxiliary strain field (for Til) within an 
C element. i=1..8. j=1..6. 

ASTRN13(i,j) Array for the auxiliary strain field (for T13) within an 

element. i=1..8. j=1..6. 

ASTRS(i,j) Array for the auxiliary stress field (for Til) within an 
element. i=1..8. j=1..6. 

ASTRS13(i,j) Array for the auxiliary stress field (for T13) within an 

element. i=1..8. j=1..6. 

AYT Total cross sectional area on the far end. 

COORD (i,j) Array for the j-th local cylindrical coordinates of the i-th 
node of the element. i=1..20. 

C00RDL(i,j) Array for the j-th local Cartesian coordinates of the i-th 

node of the element. i=1..20. 

CPL(i,j) Array for the 6x6 compliance matrix. 

C CPLR(i,j) Array for the 5x5 reduced compliance matrix. 

C CRKL 1/2 crack length. 

C DISP(i,j) Aarray for the displacements of the i-th node, i = 1..20. j = 1 . .3 
C representing UX, UY, and UZ, respectively. 

C EM ( i) Array for 3 Young's moduli. 

C ES Element thickness, or the length of "delta". 

C EZZ(i) Array for the average Ezz strains within an element. 

C EZZA(i) Array for the average Ezz strains of an element layer. 

C EX Young's modulus of the isotropic material. 

C F Area under the s-function curve. 

C FF Magnitude of the uniform force on the crack front. 

C FYT Total nodal force on the far end. 

C GM(i) Array for 3 shear moduli. 

C IRING Number of the element rings to be computed. 

C K2(i) Array for the square of the Stroh normalization factors. 

C LOAD Loading ID number. l=fixed displacement, 2=fixed stress. 

C MAT Material ID number. l=anisotropic, 2=isotropic. 

C MIE(i,j) Array for the global element number at the j-th location of the 

C i-th layer (segment) . 

C MNE(i,j,k) Array for the global node number at the j-th location of the 
C i-th layer (segment) . k=1..20 is the local node number. 

C MU ( i , j ) Array for the Storh eigenvalues. 

C N180 Number of the elements per 180 degrees. 

C N90 Number of the elements per 90 degrees. 

C NSEG Number of the element layers (segments) along crack front. 

C PL Location for the integration points. 

C PLX 1/2 panel width. 

C PLZ Panel thickness. 

C PR Poisson's ratio of the isotropic material. 

C PRM(i) Array for 3 Poisson's ratios. 

C RG ( i) , SG ( i) , TG ( i) Arrays for the natural coordinates of Gaussian 
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C integration points, i=1..8. 

C SUP The 1-1 component of the reduced compliance, all'. 

C SIFl(i, j) Array for the value of the stress intesity factor in the i-th 
C element layer of the j-th ring. 

C SIFlN(i,j) Array for the value of the normalized stress intesity factor in 
C the i-th element layer of the j-th ring. 

C SIF1NAVG ( i , 2 ) The average normalized stress intensity factor. 

C SINF Equivalent far-field stress. 

C STRN(i,j) Array for the FE result of strains within an element. The 
C i-th row stores strain components at the i-th Gaussian 

C integration point, i=1..8. Strain components from 1st to 6th 

C column are [Exx, Eyy , Ezz , Exy, Eyz , Exz ] , respectively. 

C STRS(i,j) Array for FE result of stresses within an element. See "STRN" 

C for similar definition. 

C Tll(i,j) Array for the value of the Til stress in the i-th element layer 
C of the j-th ring. 

C TllN(i,j) Array for the value of the normalized Til stress in the i-th 
element layer of the j-th ring. 

T13(i,j) Array for the value of the T13 stress in the i-th element layer 
of the j-th ring. 

T13N(i,j) Array for the value of the normalized T13 stress in the i-th 
element layer of the j-th ring. 

T33(i,j) Array for the value of the T33 stress in the i-th element layer 
of the j-th ring. 

T33N(i,j) Array for the value of the normalized T33 stress in the i-th 
element layer of the j-th ring. 

TllNAVG(i,2) The average normalized Til stress in the i-th layer. 

T13NAVG(i,2) The average normalized T13 stress in the i-th layer. 

T33NAVG(i,2) The average normalized T33 stress in the i-th layer. 

THICK(i,j) Array for Z-coordinate ( j = 1 global; j=2 normalized) of the 

center of the i-th segment, i.e. the "s" corresponding to the 
local I (s) . 

UMINU Difference of mul and mu2 (Stroh eigenvalues) . 

UPLUS Sum of mul and mu2 (Stroh eigenvalues) . 

UPROD Product of mul and mu2 (Stroh eigenvalues) . 

VEDI(i,j) Aarray for the value of local I(s) (equivalent domain 
integral of the i-th segment) in the j-th ring. 

VII (i,j) Aarray for the value of local I(s) (interaction integral 

C of the i-th segment) for Til in the j-th ring. 

C VI2(i,j) Aarray for the value of local I(s) (interaction integral 

C of the i-th segment) for T13 in the j-th ring. 

C X0J(i,j,k) Array for the value of the equivalent domain integral of 
C each element. i=l..IRING; j=l..NSEG. 

C XlJ(i,j,k) Array for the value of the interaction integral of each 
C element for Til. i=l..IRING; j=l..NSEG. 

C X2J(i,j,k) Array for the value of the interaction integral of each 
C element for T13. i=l..IRING; j=l..NSEG. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX (8) MU (3,2) , K2 , P , Q , UPLUS , UMINU, UPROD 
PARAMETER (FF = 1) 

CHARACTER HEAD*4 0 , RING* 7 
DIMENSION RING (6) 

DATA RING /'RING #2 ' , ' RING #3', 'RING #4 ' , ' RING #5', 'RING #6 ' , ' RING 
+ #7'/ 

COMMON PI , PL , HEAD 

COMMON /AUX/ ASTRS(8,6) ,ASTRN(8,6) ,ADISP(20,3) 

COMMON /AUX13/ ASTRS13 (8,6) , ASTRN13 (8,6) , ADISP13 (20,3) 

COMMON /FESOL/ DISP (2 1 , 3 ) , STRS ( 8 , 6 ) , STRN ( 8 , 6 ) 

COMMON /FEMOD/ MIE ( 50 , 24 ) , MNE ( 50 , 24 , 2 0) , COORD (21 , 3 ) , COORDL ( 2 0 , 3 ) 
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COMMON /MATL/ CPL ( 6 , 6 ) , CPLR ( 5 , 5) , EX , PR 

COMMON /MUS/ UPLUS,UMINU,UPROD,B,D,K2 (3) ,P(2,2) ,Q(2,2) 

COMMON /GIP/ RG ( 8 ) , SG ( 8 ) , TG ( 8 ) 

DIMENSION EM ( 3 ) , GM { 3 ) , PRM ( 3 ) , CPLN ( 3 ) 

DIMENSION VEDI (50 , 6 ) , VII ( 50 , 6) , X0 J (4 , 50 , 24 ) , XI J ( 4 , 50 , 24 ) 

DIMENSION EZZ (24) , EZZA(50) , THICK (50, 3) , SIF1 ( 50 , 6 ) , Til ( 50 , 6) 

DIMENSION SIF1N(50, 6) ,T11N(50, 6) 

DIMENSION VI 2 (50,6) , X2 J (4 , 50 , 24 ) ,T13 (50,6) ,T13N(50,6) 

DIMENSION T33 (50, 6) ,T33N(50, 6) 

DIMENSION SIF1NAVG (50,2) , T11NAVG (50 , 2 ) , T13NAVG (50 , 2 ) , T33NAVG (50 , 2) 

OPEN ( 1 , FILE=' sts . dat ', STATUS= ' unknown' ) ! stress data file 

OPEN (2 , FILE=' stn. dat ', STATUS= ' unknown' ) ! strain data file 

OPEN (3 , FILE=' dis . dat ', STATUS= ' unknown' ) ! displacement data file 

OPEN (4 , FILE=' node . dat ', STATUS= ' unknown' ) ! nodal data file 

OPEN (5 , FILE=' tsts . out ', STATUS= ' unknown' ) ! parameter output file 

OPEN ( 8 , FILE=' stn33 . dat ' , STATUS= ' unknown' ) ! crack-front strain data file 

OPEN ( 9 , FILE=' par . dat ', STATUS= ' unknown' ) ! parameter data file 

OPEN ( 12 , FILE= ' KIN . out ' , STATUS =' unknown' ) 

OPEN ( 13 , FILE= ' TUN . out ' , STATUS= ' unknown' ) 

OPEN (24 , FILE= ' T13N . out ' , STATUS= ' unknown' ) 

OPEN (25, FILE= ' T33N . out ' , STATUS =' unknown' ) 

PI = ACOS (-1.0) 

PL = 1 . 0/SQRT (3 . 0) ! location for integration points 

C. . .Read parameters: 

READ (9, ' (I1,TR4,I1) ' ) MAT, LOAD 
READ ( 9 , ' ( //// /TR1 1 , E16 . 9 ) ' ) AYT 

IF (MAT .EQ. 1) WRITE (5 ,'( "Material : Anisotropic")') 

IF (MAT .EQ. 2) WRITE ( 5 ,'( "Material : Isotropic")') 

IF (LOAD .EQ. 1) WRITE (5,' ("Far-field loading: Fixed Displacement 

+ ") ' ) 

IF (LOAD .EQ. 2) WRITE ( 5 , ' ( "Far - f ield loading: Fixed Load")') 

READ (9,' (/////TR11 , E16 . 9) ' ) FYT 

IF (LOAD .EQ. 1) WRITE ( 5 , ' ( A23 , IP , E16 . 9) ' ) 'Total Nodal Force FY 
+= ' , FYT 

WRITE (5, ' (A22,E14. 7) ' ) 'Cross-sectional area =',AYT 
READ (9, ' (/////TR12 , 12) ' ) NSEG 
WRITE (5, ' (A17, 13) ' ) 'No. of Segments =',NSEG 
READ (9, ' (/////TR11,F11.9) ' ) PLX 

WRITE ( 5 , ' (A18 , IP, E13 . 6 , A3 ) ' ) '1/2 Panel width = ',PLX, 'in' 

READ (9, ' (/////TR11,E16.9) ' ) PLZ 

WRITE ( 5 , ' (A23 , E13 . 6 , A3 ) ' ) 'Full Panel Thickness = ',PLZ,'in' 

GOTO (1,6) MAT ! read mat ' 1 properties based on the mat ' 1 type 

1 READ (9,3) ( EM (K) , K=1 , 3 ) , ( PRM (K) , K=1 , 3 ) , (GM (K) , K=1 , 3 ) 

3 FORMAT U////9 (TR28 , E16 . 9/) ) 

4 FORMAT (3 (E14 . 10, 5X) / ,3 (Fll . 10, 8X) / ,3 (E14 . 10, 5X) ) 

WRITE (5,5) (EM(K) ,K=1,3) , (PRM(K) ,K=1,3) , (GM(K) ,K=1,3) 

5 FORMAT ( ' Ex, Ey , Ez = ' , 3 (IP, E16 . 9 , IX) , ' ps i ' , / ' NUy z , NUxz , NUxy = ', 

+3 (OP, Fll. 9, IX) ,/'Gyz,Gxz,Gxy = ' , 3 ( IP , E16 . 9 , IX) , ' psi') 

GOTO (7,8) LOAD ! determine load type 

6 READ (9,' ( /////TR11 , E15 . 8 ) ' ) EX 

WRITE (5 , ' (A21 , E15 . 8 , A4 ) ' ) "Young's Modulus Ex = ",EX,"psi" 

READ (9, ' (/////TR12,E15.9) ' ) PR 

WRITE (5,' (A22,F6.4)') "Poisson's Ratio nuxy = " , PR 
GOTO (7,8) LOAD ! determine load type 

7 READ (9, ' (/////TR11,E16.9/) ' ) APPSTN 

WRITE (5,' (A18 , IP, E13 . 6) ' ) 'Far-field strain = ', APPSTN 
SINF = ABS (FYT/AYT) ! equivalent far-field stress 

GOTO 9 
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8 READ (9, ' (/////TR12,E15.8/) ' ) APPSTS 

WRITE (5, ' (A17,1P,E13.6,A4) ' ) 'Far-field load = ' , APPSTS psi ' 

SINF = APPSTS 

9 WRITE (5,' (A30 , IP, E13 . 6 , A4) ' ) 'Equivalent far-field stress = 

+SINF, 'psi' 

READ (9, ' (////TR12 , E13 . 6) ' ) CRKL 

WRITE (5, ' (A18,1P,E13.6,A3) ' ) '1/2 crack length = ', CRKL, 'in' 

READ (9, ' (/////TR12 , II) ' ) N90 
N180 = N90*2 

ARC = 180/N180 ! arc length (in degree) of a single element 

WRITE (5,' (A30,I2)') 'Elements in 180-degree span = ' ,N180 
WRITE (5,' (A28 , F5 . 1 , A8) ' ) 'Arc length of each element =',ARC, 

+ ' Degrees' 

READ (9 , ' ( /// //TR12 , II) ' ) IRING 

WRITE (5,' (A33,I1)') 'Number of rings to be computed = ', IRING 
GOTO (11,36) MAT 

C... Evaluate the full and reduced compliance matrices for the anisotropic 
C material : 

11 CALL STRUC ( EM , GM , PRM) 

WRITE (5,12) 'Structural Compliance Matrix [S] = ', 

+( ( CPL ( I , J) , J=1 , 6) ,1=1,6) 

12 FORMAT ( / A50/ ,6(6(1P,E13.6, IX) / ) ) 

WRITE (5,13) 'Reduced Structural Compliance Matrix "S_0" = ', 

+( (CPLR(I, J) , J=l,5) , 1=1,5) 

13 FORMAT ( / A46/ ,5(5(1P,E13.6, IX) / ) ) 

WRITE ( 5 , ' (A18 , IP , E12 . 5 ) ' ) " Sqrt (s22 ' /sll ' ) = ", 

+SQRT ( CPLR (2,2) /CPLR(1,1) ) 

C... Solve for Storh eigenvalues: 

CALL STROH(MU) 

WRITE (5,21) ' Stroh Eigenvalues =',( (MU(M,N) ,N=1, 2) ,M=1, 3) 

21 FORMAT ( / A2 0 / , 3 ( 2 ( ' (',1P,E13.6,' + ' , IP , E13 . 6 , ' i )',5X)/)) 

WRITE (5,22) 'p: ' , (M, (P (M, N) , N=1 , 2 ) , M=1 , 2) 

WRITE (5,22) 'q: ' , (M, (Q(M,N) ,N=1,2) ,M=1,2) 

22 FORMAT (A3/ , 2 ( II , 2X, 2 ( ' ( ' , IP, E13 . 6 , ' + ' , IP , E13 . 6 , ' i ) ' , 5X) / ) ) 

WRITE (5,23) 'k A 2: ' , (M, K2 (M) ,M=1, 3) 

23 FORMAT (A5/, 3 (I1,2X, ' ( ' , IP, E13 .6 , ' + ' , IP , E13 . 6 , ' i )'/)) 

C ... Calculate L_22 value, for later calculation of K_I : 

UPLUS = MU (1,1) +MU (2,1) 1 sum of 2 eigenvalues 

UPROD = MU (1 , 1) *MU(2 , 1) ! product of 2 eigenvalues 

A1 = DREAL (UPLUS) ! real part of UPLUS 

A2 = DREAL (UPROD) ! real part of UPROD 

B = DIMAG (UPLUS) ! imaginary part of UPLUS 

D = DIMAG (UPROD) ! imaginary part of UPROD 

AB = A1*D-A2*B 

BL22I = CPLR (1,1) *AB ! L_inverse_22 
CPLN(l) = CPL (1,1) I compliance sll 

CPLN ( 2 ) = CPL (1,3) 1 compliance sl3 

CPLN ( 3 ) = CPL (3,3) ! compliance s33 

BL2 1 = -D/ (CPLR (1,1)* (B*AB-D* *2 ) ) ! L2 1 

BL22 = B/ ( CPLR (1,1)* (B*AB-D* *2 ) ) ! L22 

CPLN ( 1) = CPL (1,1) I compliance sll 

CPLN (2) = CPL (1,3) ! compliance sl3 

CPLN ( 3 ) = CPL (3,3) ! compliance s33 

GOTO 40 

C...For isotropic materials: 

36 CPLN ( 1) = 1/EX 

CPLN (2) = -PR/EX 
CPLN (3) = 1/EX 
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BL22 = EX/ (2* (1-PR**2) ) 

C...Set Gauss integration point locations within an element: 

4 0 CALL GUASS 

SUP = CPLN ( 1) -CPLN (2 ) * *2/ CPLN ( 3 ) ! sll' 

C. .. Calculate the equivalent domain integral and the interaction integral in 

C every element: 

DO 300 KR = 1,IRING ! KR : counter for the rings 

DO 200 I = 1,NSEG ! I: counter for the element layers (segments) 

C... Array initialization: 

IF (KR .EQ. 1) THEN ! initialization for later average calculation 
SIF1NAVG (1,1) = 0 
T1 1NAVG (1,1) = 0 
T13NAVG (1,1) = 0 
T3 3NAVG (1,1) = 0 
END IF 

VEDI ( I , KR) = 0 
VII ( I , KR) = 0 
VI2 ( I , KR) = 0 
IF (KR .EQ. 1) THEN 
EZZA(I) = 0 
END IF 

DO 100 J = 1,N180 ! J: counter for the elements in a segment 

CALL FERST ( I , J) 

C...Set up auxiliary fields: 

61 CALL AUXT13 (MU, ARC, MAT) ! auxiliary field for T13 

GOTO (62,63) MAT 

62 CALL AUXAN (MU, ARC) ! auxiliary field for Til - anisotropic 

GOTO 70 

63 CALL AUXI (ARC) ! auxiliary field for Til - isotropic 

C. . . Evaluate equivalent domain integral and interaction integral: 

70 CALL INTEG ( KR , I , J, X0 J, XI J, X2 J) 

VEDI ( I , KR) = VEDI ( I , KR) +X0J (KR, I , J) ! equivalent domain integral 
VII (I, KR) = VII ( I , KR) +X1 J (KR, I , J) ! interaction integral for Til 

VI2 ( I , KR) = VI2 ( I , KR) +X2 J (KR , I , J) ! interaction integral for T13 

C ... Calculate average Ezz strain in every element on the crack front: 

IF (KR .EQ. 1) THEN 
CALL STRN33 (EZZ, J) 

EZZA(I) = EZZA(I) +EZZ (J) 

END IF 

100 CONTINUE 

IF (KR .EQ. 1) THEN 

EZZA(I) = EZZA(I)/N180 ! average over a layer ( segment ) 

ENDIF 

C ... Calculate T-stress and K1 using interaction- integral approach. 

ES = COORDL (5,3) - COORDL (1,3) ! element thickness 

F = ( 2 . / 3 . ) *ES 

THICK (1,1) = ( COORDL (1,3) +COORDL (5,3) ) / 2 ! local "s" coordinate 

THICK (1,2) = THICK (1,1) / ( PLZ/2 ) ! normalized thickness 

THICK (I, 3) = COORDL (1,3)/ (PLZ/2) ! local "s" coordinate at 

* bottom of the element 

VEDI ( I , KR) = 2 *VEDI ( I , KR) /F 
VII ( I , KR) = 2* VI I ( I , KR) /F 
VI2 ( I , KR) = 2 *VI2 ( I , KR) /F 
GOTO (111,112) MAT 
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C...For anisotropic materials: 

111 SIF1 ( I , KR) = SQRT (2*VEDI ( I , KR) /BL22I) 

T11(I,KR) = 1* (VII (I,KR) -EZZA(I) *CPL(1,3) /CPL(3,3) ) /CPLR(1,1) 

T13(I,KR) = VI2 (I,KR) /CPLR(5,5) 

T33 ( I , KR) = (EZZA ( I ) - (CPL (1,3) *T11 ( I , KR) +CPL (3,5) *T13 (I , KR) ) ) / 

+ CPL (3, 3) 

GOTO 113 

C...For isotropic materials: 

112 SIF1 ( I , KR) = SQRT (VEDI (I, KR) *EX/ (1-PR** * * 2) ) 

Til (I , KR) = 1* (EX/ ( ( 1 - PR* *2 ) *FF) ) * (VII ( I , KR) +FF*PR*EZZA ( I ) ) 

T13 ( I , KR) = (EX/ (2* (1+PR) ) ) *VI2 ( I , KR) 

T33 ( I , KR) = EX* EZZA ( I ) + PR* Til ( I , KR) 

113 SIF1N (I , KR) = SIF1 (I, KR) / (SINF*SQRT (PI*CRKL) ) 

TUN ( I , KR) = Til ( I , KR) /SINF 

T13N ( I , KR) = T13 ( I , KR) / SINF 
T33N ( I , KR) = T33 (I,KR) /SINF 
SIF1NAVG (1,1) = SIF1NAVG (1,1) +SIF1N ( I , KR) 

T11NAVG (1,1) = T1 1NAVG (1,1) +T11N ( I , KR) 

T13NAVG (1,1) = T13NAVG (1,1) +T13N ( I , KR) 

T33NAVG (1,1) = T33NAVG ( 1 , 1 ) +T3 3N ( I , KR) 

IF (KR .EQ. IRING) THEN ! average over rings 

SIF1NAVG (1,2) = SIF1NAVG(I, 1) /IRING 
T1 1NAVG (1,2) = T11NAVG(I, 1) /IRING 
T13NAVG (1,2) = T13NAVG(I, 1) /IRING 
T3 3NAVG (1,2) = T33NAVG (1 , 1) /IRING 
END IF 

200 CONTINUE 

300 CONTINUE 

C... Output calculated T- stresses and K1 data: 

WRITE (12,301) ' SEGMENT' , 'THICKNESS' , (RING ( K) ,K=1, IRING) , 'AVERAGE' 

WRITE (13,301) ' SEGMENT' , 'THICKNESS' , (RING ( K) ,K=1, IRING) , 'AVERAGE' 

WRITE (24,301) ' SEGMENT' , 'THICKNESS' , (RING(K) ,K=1, IRING) , 'AVERAGE' 

WRITE (25,301) ' SEGMENT' , 'THICKNESS' , (RING(K) ,K=1, IRING) , 'AVERAGE' 

301 FORMAT (A7 , IX, A9 , 4X , 4 (4X, A7 , 4X) ) 

302 FORMAT ( 14 , 4X, F9 . 6 , 6 (2X, IP, E13 . 6) ) 

DO 310 I = 1 , NSEG 

WRITE (12,302) I, THICK (1,2) , (SIF1N(I,K) , K=l, IRING) , SIF1NAVG (1,2) 

WRITE (13,302) I , THICK ( 1 , 2 ) , (TUN ( I , K) , K=1 , IRING) , T11NAVG ( 1 , 2 ) 

WRITE (24,302) I , THICK (I , 2 ) , (T13N ( I , K) , K=1 , IRING) , T13NAVG ( I , 2 ) 

WRITE (25,302) I , THICK (I , 2 ) , (T33N ( I , K) , K=1 , IRING) , T33NAVG (I , 2 ) 

310 CONTINUE 

999 STOP 
END 

SUBROUTINE STRUC (YOUNG, SHEAR, PSNR) 

*********★*★★*****★★★★★★★★*****★**★**★** 

* The subroutine computes the full and reduced compliance matrices. * 

•k k 

* Subroutine input: YOUNG, SHEAR, PSNR * 

* Subroutine output: CPL,CPLR * 

k k 

* Definition of local variables: * 

* C(i,j) Array for the stiffness matrix. * 

* C0(i,j) Array for the reduced stiffness matrix. * 

* PSNR(i) Array for 3 Poisson's ratios. * 

* SHEAR (i) Array for 3 shear moduli. * 
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* YOUNG (i) Array for 3 Young's moduli. 

★ ★★★★★★★★★★★★★★★★★★★★★it 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
COMMON /MATL/ CPL ( 6 , 6 ) , CPLR ( 5 , 5) , EX, PR 
DIMENSION YOUNG (3) , SHEAR (3) ,PSNR(3) 
DIMENSION C(6, 6) , CO (5, 5) 


★ 

★ 


C. . .Construct the 6x6 compliance matrix [S] : 
DO 915 I = 1,3 
CPL (1,1) = 1/YOUNG (I) 

CPL ( 1+3 , 1+3 ) = 1 / SHEAR ( I ) 


915 


CONTINUE 
CPL (2,1) 
CPL (1 , 
CPL (3 , 
CPL ( 1 , 
CPL (3, 
CPL (2, 


, 2 ) 

,D 

,3) 

, 2 ) 

,3) 


- (PSNR ( 3 ) /YOUNG ( 1 ) ( 
CPL (2 ,1) 

- (PSNR (2) /YOUNG (1) ( 
CPL (3 , 1) 

- ( PSNR ( 1 ) /YOUNG (2) ( 
CPL (3,2) 


C... Evaluate the 6x6 stiffness matrix [C] , the 5x5 reduced stiffness and the 
C compliance matrices: 

CALL DLINRG ( 6 , CPL, 6 , C, 6 ) ! IMSL Library for matrix inversion 

DO 920 I = 1,5 
DO 920 J = 1,5 

IF ((I .NE. 3) .AND. (J .NE. 3)) THEN 
CO (I, J) = C ( I , J) 

ELSEIF ( (J .EQ. 3) .AND. (I .NE. 3)) THEN 
CO (I, J) = C(I,6) 

ELSEIF ((I .EQ. 3) .AND. (J .NE. 3)) THEN 
CO (I, J) = C (6 , J) 

END IF 

920 CONTINUE 

CO (3,3) = C ( 6 , 6) 

CALL DLINRG (5 , CO , 5 , CPLR, 5) ! IMSL Library for matrix inversion 

RETURN 

END 


SUBROUTINE STROH(U) 

★ ★★★★★*★★★★★★★★★★★★★★★★★★★★★★*★*★★★★★★★★ 

* The subroutine determines the Stroh eigenvalues from the reduced * 

* structural compliance matrix. * 

* ★ 

* Subroutine input: CPLR * 

* Subroutine output: U,K2,P,Q * 

* ★ 


* Definition of local variables: * 

* Al(i) Array for the 5 coefficients in solving the in-plane * 

* characteristic equation. * 

* U(m,n) Array for the Stroh eigenvalues. m=1..3, 3 pairs of eigenvalues* 

* and m=3 indicates the out-of-plane eigenvalue. n=2 is the * 

* conjugate of n=l. * 

* Z(i) Array for the 4 roots of the in-plane characteristic equation. * 

* ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX ( 8 ) U,Z,D2,K2,P,Q, UPLUS , UMINU, UPROD 
COMMON /MATL/ CPL ( 6 , 6 ) , CPLR ( 5 , 5) , EX , PR 

COMMON /MUS/ UPLUS, UMINU, UPROD, B,D,K2 (3) ,P(2,2) ,Q(2,2) 
DIMENSION U ( 3 , 2 ) , A1 ( 5 ) , Z ( 4 ) 


C... Solve for 4 in-plane eigenvalues: 
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c. 

c 


A1 (1) 

= CPLR (2, 2) 

! coefficient 

for 

mu A 0 : "s'_22" 

A1 (2) 

= -2 *CPLR (2,3) 

! coefficient 

for 

mu A l : " - 2 * s ' _2 6 " 

A1 (3) 

= 2*CPLR (1,2) +CPLR(3, 3) 

! coefficient for mu A 2 

A1 (4) 

= -2 *CPLR (1,3) 

! coefficient 

for 

mu A 3 : " -2 *s ' _16 " 

A1 (5) 

= CPLR (1,1) 

! coefficient 

for 

mu A 4 : "s' 11" 


..Use an IMSL Library to find the roots of a polynomial with real 
coefficients by Laguerre's method: 


CALL DZPLRC ( 4 , A1 , Z) 




U(l,l) = Z(l) 

! "mul " 



U ( 1 , 2 ) = Z (2) 

I "mul_bar " 



U(2 , 1) = Z (3) 

! "mu2 " 



U(2 ,2) = Z (4) 

! "mu2_bar " 



UMINU = U(l, 1) -U (2 , 1) 

! "mul-mu2" 



lve for 2 out -of -plane 

eigenvalues : 



A2 = CPLR (5, 5) 

! coefficient 

for 

A ~ 

mu 2 

B2 = -2 *CPLR (4,5) 

! coefficient 

for 

mu 1 

C2 = CPLR (4, 4) 

! coefficient 

for 

o 

< 

3 

B 

D2 = B2 * * 2 - 4 * A2 * C2 





U ( 3 , 1 ) = (-B2+SQRT (D2) ) / (2*A2) ! "mu3" 

U ( 3 , 2 ) = ( -B2 - SQRT (D2 ) ) / (2 *A2 ) ! "mu3_bar" 

C ... Calculate p, q, and normalization factors k: 


" s ' _5 5 " 

" - 2 * s ' _4 5 " 
"s' 44" 


DO 935 N = 1,2 

P (N, 1) = CPLR (1,1) *U (N, 1) **2-CPLR (1,3) *U (N, 1) +CPLR (1,2) 
P (N, 2 ) = CONJG ( P (N , 1 ) ) 

Q (N, 1) = CPLR (1,2) *U (N , 1 ) - CPLR {2,3) +CPLR (2,2)/U(N,l) 

Q (N, 2 ) = CON J G ( Q ( N , 1 ) ) 

K2 (N) = 1/ (2* (Q(N, 1) -U(N,1) *P(N, 1) ) ) 

935 CONTINUE 

K2 ( 3 ) = 1/ (2* (CPLR(4,4) /U(3, 1) -CPLR(4,5) ) ) 

990 RETURN 
END 


! "pi" 

! "pl_bar " 

! " ql " 

! "ql_bar" 

! kl A 2 , k2 a 2 

! k3 A 2 


SUBROUTINE GUASS 

**************************************** 

* The subroutine defines the Gaussian integration point locations using * 

* reduced 8-point rule. * 

•k * 

* Subroutine input: none * 

* Subroutine output: RG,SG,TG * 

* * 

* Definition of local variables: * 

* PL1 Integration point location. * 

* RI,SI,TI Coeffecients for the integration point locations. * 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /GIP/ RG(8) ,SG(8) ,TG(8) 

DIMENSION SI (8) ,TI (8) ,RI (8) 

DATA RI /-l, 1, 1, -1, -1,1,1, -1/ 

DATA SI /-l, -1,1, 1,-1, -1,1,1/ 

DATA TI /-l, -1,-1, -1,1, 1,1,1/ 


PL1 = 1/SQRT(3.0) ! point location for integration 

C...Set up 2x2x2 reduced Gaussian integration points: 

DO 10 K = 1,8 
RG (K) = PL1*RI(K) 

SG (K) = PL1*SI(K) 

TG (K) = PL1*TI(K) 

10 CONTINUE 
RETURN 
END 



95 


SUBROUTINE FERST (KS , KE) 

★ ★★★★★★★★★★★★★★★★★★★★★★★★it************** 

* The subroutine reads strain, stress, and nodal data from ANSYS results.* 

* Those data are stored in different arrays for later uses. * 

★ ★ 

* Subroutine input: data files 1,2, 3, 4 * 

* Subroutine output: DISP, STRS, STRN, MIE,MNE , COORD, COORDL * 

* ★ 

* Definition of local variables: * 

* KE The element number in the layer (segment) . * 

* KS The element layer (segment ) number along the crack front. * 

* THETA Angular coordinate (in radians) of the nodal points. * 

* ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

CHARACTER HEAD*40 
COMMON PI , PL , HEAD 

COMMON /FESOL/ DISP (2 1 , 3 ) , STRS ( 8 , 6 ) , STRN ( 8 , 6 ) 

COMMON /FEMOD/ MIE (50,24) , MNE ( 50 , 24 , 2 0) , COORD (2 1,3) , COORDL ( 2 0 , 3 ) 

C...Read element -nodal data and sort by local node numbers: 

READ (4,1001) MIE (KS, KE) , (MNE (KS , KE , K) , K=1 , 20) 

1001 FORMAT (///////////////////2X,I6, TR2 1 , 8 16/TR29 , 8 16/TR2 9 , 

+416//////////////////) ! ANSYS 5.5 format 

II = 0 

1002 READ (4,1003) NK, (COORD (21 , J) , J= 1 , 3 ) ! NK: node no. in the data file 

1003 FORMAT ( IX, 17 , 2X, 3E12 . 5) 

II = II+l 

DO 1010 I = 1,20 ! I: local node no. 

NI = MNE ( KS , KE , I ) ! NI : global node no. 

IF (NI .EQ. NK) THEN 
DO 1005 J =1,3 
COORD ( I, J) = COORD (21, J) 

1005 CONTINUE 

IF (II .LT. 20) GOTO 1002 
END IF 

1010 CONTINUE 


C... Convert local cylindrical nodal coordinates to local Cartesian coordinates: 
DO 1020 I = 1,20 ! I: local node no. 

THETA = COORD (1,2) *PI/ 180 

! xl -coordinate 
! x2 -coordinate 
! x3 -coordinate 



COORDL (1,1) = 

COORD ( I, 1) 

*COS (THETA) 


COORDL (1,2) = 

COORD ( I, 1) 

* SIN (THETA) 


COORDL (1,3) = 

COORD (1,3) 


1020 

CONTINUE 



C. . . 

Read stress data 

within an 

element : 


READ (1,1021) 

( STRS ( 1 , J) , 

J=1 , 6) 

1021 

FORMAT (//// 9X 

, 6E12 . 5 ) 

I ANSYS 

1023 

FORMAT ( 9X , 6E12 . 5) 



DO 1030 1=2, 

8 



READ (1,1023) 

( STRS ( I , J) 

, J=l,6) 

1030 

CONTINUE 



C. . . 

Read strain data 

within an 

element : 


READ (2,1021) 

( STRN ( 1 , J) , 

J=1 , 6) 


DO 1040 1=2, 

8 



READ (2,1023) 

( STRN ( I , J) 

, J=l, 6) 

1040 

CONTINUE 




C...Sort stress and strain data in consistent with the stress and strain 
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C vectors, i.e. [sigma] = [Sxx, Syy, Szz , Syz , Sxz , Sxy] and similar to strains: 
DO 1050 I =1,8 
SXY = STRS (1,4) 

STRS (1,4) = STRS (I, 5) 

STRS (I, 5) = STRS (I, 6) 

STRS (1,6) = SXY 
EXY = STRN (1,4) 

STRN (1,4) = STRN (I, 5) 

STRN (I, 5) = STRN (I, 6) 

STRN (I, 6) = EXY 
1050 CONTINUE 


C...Read displacement data within an element and sort by local node numbers: 

1051 READ ( 3 , ' ( A4 0 ) ' ) HEAD 
II = 0 

1052 READ (3,1053) NK, (DISP (21 , J) , J=1 , 3 ) ! NK: node no. in the data file 

1053 FORMAT ( IX, 17 , IX , 3E12 . 5 ) 

II = II+l 

DO 1060 I =1,20 | I: local node no. 

NI = MNE ( KS , KE , I ) ! NI : global node no. 

IF (NK .EQ. NI) THEN 
DO 1055 J =1,3 
DISP ( I , J) = DISP (2 1 , J) 

1055 CONTINUE 

IF (II .LE. 19) GOTO 1052 
END IF 

1060 CONTINUE 
1100 RETURN 
END 

SUBROUTINE AUXT13 (U, AL, MAT) 

**************************************** 

* The subroutine calculates the auxiliary stress, strain, and displacement* 

* fields under a line load f3 applying on the crack front, in order to * 

* calculate the T13 stresses later. * 


* Subroutine input: U, AL , MAT 

* --- Subroutine output: ASTRS13 , ASTRN13 , ADISP13 

* 


★ 

★ 

★ 

★ 

■k 

k 

k 

k 

k 


Definition of local variables: 


AL 

R 

RE 

SR3 (i) 
THE TAD 

U(i) 


Arc length (in degrees) of a single element. 

Radial coordinate of the integration or nodal points. 

Element size in r-direction. 

r-3 stresses on the 8 integration points. 

Angular coordinate (in degrees) of the integration or nodal 
points . 

Stroh eigenvalues, i=l,2,3. 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX (8) U ( 3 , 2 ) , P , Q , XI3 , OMEGA3 , K2 , UPLUS , UMINU, UPROD 
PARAMETER (FF = 1) 

COMMON PI , PL , HEAD 

COMMON / AUX1 3 / ASTRS13 (8,6) , ASTRN13 (8,6) , ADISP13 (20,3) 

COMMON /FEMOD/ MIE ( 50 , 24 ) , MNE ( 50 , 24 , 2 0) , COORD (21 , 3 ) , COORDL (2 0 , 3 ) 
COMMON /MATL/ CPL ( 6 , 6 ) , CPLR ( 5 , 5) , EX, PR 

COMMON /MUS/ UPLUS, UMINU, UPROD, B,D,K2 (3) ,P(2,2) ,Q(2,2) 

DIMENSION JR ( 8 ) , JS ( 8) 

DIMENSION SR3 ( 8) 

DOUBLE PRECISION MU 

DATA JR /l, 1,-1, -1,1,1, -1,-1/ 


k 


k 
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DATA JS /- 1, 1, 1, -1, -1,1,1, -1/ 

RE = COORD (1,1) - COORD (4,1) 

GOTO (1302,1351) MAT ! determine anisotropic or isotropic material 
C...For anisotropic material: 

1302 MU = (CPLR (4,4) *CPLR (5,5) -CPLR (4,5) **2) ** (-0.5) 

DO 1320 L = 1,8 ! L: counter for the integration points 

R = COORD (L, 1) - (1-PL) * (RE/2) * JR (L) 

THETAD = COORD(L,2) - (1-PL) * (AL/2) *JS (L) 

XI3 = COSD (THETAD) +U (3,1) *SIND (THETAD) 

0MEGA3 = (-SIND (THETAD) +U(3 , 1) *COSD (THETAD) ) /XI3 


A1 = REAL ( K2 ( 3 ) *OMEGA3 ) 

SR3 (L) = FF*A1/ (PI*MU*R) ! "Sr3" 

C... Compute Cartesian stress components: 

ASTRS13 (L, 1) = 0 ! "Sxx" 

ASTRS13 ( L , 2 ) = 0 ! "Syy" 

ASTRS13 ( L , 3 ) = 0 ! "Szz" 

ASTRS13 ( L , 4 ) = SIND (THETAD) *SR3 (L) ! "Syz" 

ASTRS13 ( L , 5 ) = COSD (THETAD) *SR3 (L) ! "Sxz" 

ASTRS13 ( L , 6 ) = 0 ! "Sxy" 

C... Compute engineering strains: 

DO 1310 I = 1,6 
ASTRN13 (L, I) = 0 
DO 1310 J = 1,6 


ASTRN13 (L , I ) = ASTRN13 (L , I) +CPL ( I , J) *ASTRS13 (L, J) 

1310 CONTINUE 

C... Convert engineering strains to tensor ial strains: 

ASTRN13 (L, 4 ) = ASTRN13 (L , 4 ) /2 
ASTRN13 (L, 5) = ASTRN13 (L , 5 ) /2 
ASTRN13 (L, 6) = ASTRN13 (L, 6) /2 
1320 CONTINUE 
C... Compute displacements: 

DO 1330 I = 1,20 ! I: counter for the local nodes 

R = COORD (1,1) 

THETAD = COORD (I, 2) 

XI3 = COSD ( THETAD) +U (3,1) *SIND ( THETAD) 

ADISP13 (1,1) = 0 ! "ul" 

ADISP13 (1,2) = 0 ! "u2 " 

ADISP13 (1,3) = -FF* ( LOG (R) +REAL (LOG (XI 3 ) ) ) / ( 2 *PI*MU) !"u3" 

1330 CONTINUE 
GOTO 1400 

C...For isotropic material: 

1351 DO 1370 L = 1,8 

R = COORD (L, 1) - (1-PL) * (RE/2) *JR(L) 

THETAD = COORD(L,2) - (1-PL) * (AL/2) * JS (L) 

SR3 (L) = -FF/(2*PI*R) ! "Sr3" 

C... Compute Cartesian stress components: 


ASTRS13 (L, 1) = 0 1 "Sxx" 

ASTRS13 (L , 2 ) = 0 ! "Syy" 

ASTRS13 (L , 3 ) = 0 ! "Szz" 

ASTRS13 (L, 4 ) = -FF*SIND (THETAD) / (2 *PI*R) ! "Syz" 

ASTRS13 (L , 5 ) = -FF*COSD (THETAD) / (2 *PI*R) ! "Sxz" 

ASTRS13 ( L , 6 ) = 0 ! "Sxy" 

C... Compute tensorial strains: 


ASTRN13 (L , 1 ) = (ASTRS13 (L, 1) -PR* (ASTRS13 (L, 2) +ASTRS13 (L, 3) ) ) /EX ! "Exx" 
ASTRN13 (L,2) = (ASTRS13 ( L, 2 ) - PR* (ASTRS13 (L , 1) +ASTRS13 (L , 3 ) ) ) /EX ! "Eyy" 
ASTRN13 ( L , 3 ) = (ASTRS13 (L, 3 ) - PR* (ASTRS13 (L , 1) +ASTRS13 (L , 2 ) ) ) /EX ! "Ezz" 
ASTRN13 ( L , 4 ) = ( 1+PR) /EX*ASTRS13 (L , 4 ) ! "Eyz" 
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ASTRN13(L,5) = ( 1+PR) /EX*ASTRS13 (L , 5 ) 
ASTRN1 3 ( L , 6 ) = ( 1+PR) /EX*ASTRS13 (L , 6 ) 
1370 CONTINUE 

DO 1380 I = 1,20 
R = COORD (1,1) 

ADISP13 (1,1) = 0 
ADISP13 (1,2) = 0 

ADISP13 (1,3) = -FF* (1+PR) *LOG(R) / (PI*EX) 
1380 CONTINUE 


1 "Exz" 
1 "Exy " 


! "ul" 
! "u2 " 
1 "u3 " 


1400 RETURN 
END 


SUBROUTINE AUXAN(U,AL) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* The subroutine calculates the auxiliary stress, strain, and displacement* 

* fields under a line load fl applying on the crack front, for an anisotropic* 

* material according to the Stroh formalism with the normalization factor k. * 

* The auxiliary fields are later used to determine Til stresses. * 

* * 

* Subroutine input: U,AL * 

* Subroutine output: ASTRS , ASTRN, ADISP * 

* * 


•k 

k 

k 

k 

k 

k 

k 

k 
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Definition of local variables: 

AL Arc length (in degrees) of a single element. 

R Radial coordinate of the integration or nodal points. 

RE Element size in r-direction. 

SRR(i) r-r stresses on the 8 integration points. 

THETAD Angular coordinate (in degrees) of the integration or nodal 

points . 

U(i) Stroh eigenvalues, i=l,2,3. 

************************************* 


* 


* 

k 

k 

k 

k 

k 

k 

k 

k 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX (8) U ( 3 , 2 ) , P , Q , XI ( 2 ) , OMEGA (2) , K2 , UPLUS , UMINU, 

+UPR0D , DL (2,2) ,S(2,2) 

PARAMETER (FF = 1) 

COMMON PI, PL, HEAD 

COMMON /AUX/ ASTRS (8,6) , ASTRN (8,6) , ADISP (20,3) 

COMMON /FEMOD/ MIE ( 50 , 24 ) , MNE ( 50 , 24 , 2 0) , COORD ( 2 1 , 3 ) , COORDL ( 2 0 , 3 ) 

COMMON /MATL/ CPL ( 6 , 6 ) , CPLR ( 5 , 5) , EX, PR 

COMMON /MUS/ UPLUS, UMINU, UPROD,B,D,K2 (3) ,P(2,2) , Q (2 , 2 ) 

DIMENSION JR ( 8 ) , JS ( 8) 

DIMENSION SRR ( 8 ) 

DATA JR /l, 1,-1, -1,1,1, -1,-1/ 

DATA JS /-1,1, 1,-1, -1,1, 1,-1/ 


RE = COORD (1,1)- COORD (4,1) 

DO 1320 L = 1,8 ! L: counter for the integration points 

R = COORD (L, 1) - (1-PL) * (RE/2) *JR(L) 

THETAD = COORD ( L, 2 ) - ( 1-PL) * (AL/ 2 ) * JS (L) 

XI (1) = COSD ( THETAD) +U( 1,1) *SIND (THETAD) 

XI (2) = COSD (THETAD) +U (2,1) *SIND (THETAD) 

OMEGA ( 1 ) = (-SIND (THETAD) +U ( 1 , 1) *COSD (THETAD) ) /XI (1) 

OMEGA ( 2 ) = (-SIND (THETAD) +U(2 , 1) *COSD (THETAD) ) /XI (2) 

DL (1,1) = K2 (1) *U ( 1 , 1) **2*OMEGA ( 1 ) +K2 (2 ) *U (2 , 1) **2 *OMEGA (2 ) 

DL (1,2) = -K2 (1) *U(1, 1) *OMEGA(l) -K2 (2) *U(2 , 1) *OMEGA(2) 

DL (2,1) = DL (1,2) 

DL (2,2) = K2 ( 1) *OMEGA ( 1 ) +K2 (2 ) * OMEGA (2 ) 

A1 = REAL (COSD (THETAD) * (B*DL(1, 1) +D*DL (1,2) ) +SIND (THETAD) * 

+ (B*DL (2,1) +D*DL (2,2))) 
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SRR(L) = FF*CPLR(1 , 1) *A1/ (PI*R) ! "Srr" 

C... Compute Cartesian stress components: 

ASTRS(L,1) = FF*CPLR (1,1) *COSD ( THETAD) **2/ (PI *R) *A1 ! "Sxx" 

ASTRS(L,2) = FF*CPLR (1,1) *SIND ( THETAD) **2/ (PI *R) *A1 ! "Syy" 

ASTRS (L / 3 ) = -FF*CPLR (1,1) * (CPL (1,3) *COSD (THETAD) **2+CPL (2,3) 

+ *SIND (THETAD) **2 ) / (PI*R*CPL (3,3) ) *A1 ! "Szz" 

ASTRS (L, 4) = 0 1 "Syz" 

ASTRS (L, 5) = 0 1 "Sxz" 

ASTRS ( L, 6 ) = FF*CPLR (1,1) *SIND ( THETAD) *COSD (THETAD) / (PI*R) *A1 ! "Sxy" 

C... Compute engineering strains: 

DO 1310 I = 1,6 
ASTRN (L, I) = 0 
DO 1310 J = 1,6 

ASTRN (L, I) = ASTRN (L , I ) +CPL ( I , J) *ASTRS (L , J) 

1310 CONTINUE 

C... Convert engineering strains to tensorial strains: 

ASTRN (L, 4 ) = ASTRN (L , 4 ) /2 
ASTRN (L, 5) = ASTRN (L, 5) /2 
ASTRN (L, 6) = ASTRN (L , 6 ) /2 
1320 CONTINUE 
C... Compute displacements: 

DO 1330 I =1,20 ! I: counter for the local nodes 

R = COORD (1,1) 

THETAD = COORD (I, 2) 

XI (1) = COSD ( THETAD) +U(1, 1) *SIND( THETAD) 

XI (2) = COSD (THETAD) +U (2,1) *SIND (THETAD) 

S(l,l) = -K2 ( 1) *U (1,1)*P(1,1) *LOG ( XI ( 1 ) ) -K2 (2) *U(2, 1) *P(2, 1) * 

+ LOG (XI (2) ) 

S ( 1 , 2 ) = K2 ( 1 ) *P ( 1 , 1 ) * LOG (XI ( 1 ) ) +K2 (2 ) *P (2 , 1 ) *LOG (XI (2 ) ) 

S ( 2 , 1 ) = -K2 (1) *U(1,1) *Q(1,1) *LOG(XI(l) ) -K2 (2) *U(2,1) *Q(2,1) * 

+ LOG (XI (2) ) 

S ( 2 , 2 ) = K2 ( 1 ) *Q ( 1 , 1 ) *LOG (XI ( 1) )+K2 (2) *Q(2, 1) *LOG (XI (2 ) ) 

ADISP (1,1) = -FF*CPLR (1,1) * (B*LOG ( R) +2 *REAL (B * S ( 1 , 1 ) +D*S (1,2) )) / 

+ (2 *PI ) 1 "ul" 

ADISP (I, 2) = -FF*CPLR (1,1) * (D*LOG (R) +2*REAL(B*S (2,1) +D*S (2,2)))/ 

+ ( 2 *PI ) ! "u2 " 

ADISP (I, 3) = 0 ! "u3 " 

1330 CONTINUE 
1400 RETURN 
END 


★ 
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SUBROUTINE AUXI (AL) 

The subroutine calculates the auxiliary stress, strain, and displacement* 
fields under a line load applying on the crack front, for an isotropic * 

material. The auxiliary fields are later used to determine Til stresses. * 

* 

Subroutine input: AL * 

Subroutine output: ASTRS, ASTRN, ADISP * 

* 

Definition of local variables: * 

AL Arc length (in degrees) of a single element. * 

RE Element size in r-direction. * 

THETA Angular coordinate (in radians) of the integration or nodal * 

points. * 

THETAD Angular coordinate (in degrees) of the integration or nodal * 

points. * 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
PARAMETER (FF = 1) 
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COMMON PI , PL , HEAD 

COMMON /AUX/ ASTRS(8,6) ,ASTRN(8,6) ,ADISP(20,3) 

COMMON /FEMOD/ MIE ( 50 , 24 ) , MNE ( 50 , 24 , 2 0) , COORD (21 , 3 ) , COORDL ( 2 0 , 3 ) 
COMMON /MATL/ CPL ( 6 , 6 ) , CPLR ( 5 , 5) , EX, PR 
DIMENSION JX1 ( 8) , JX2 (8) 

DATA JX1 /l, 1,-1, -1,1,1, -1,-1/ 

DATA JX2 / - 1 , 1 , 1,-1, -1,1, 1,-1/ 

RE = COORD (1 , 1) -COORD (4 , 1) ! element size in r-direction 

DO 1310 L = 1,8 ! L: counter for the integration points 


R = COORD (L, 1) - (1-PL) * (RE/2) *JX1 (L) 

THETAD = COORD (L, 2) - (1-PL) * (AL/2) *JX2 (L) 

ASTRS (L, 1 ) = -FF*COSD (THETAD) **3/ (PI*R) ! "Sxx" 

ASTRS ( L, 2 ) = -FF*COSD (THETAD) *SIND (THETAD) **2/ (PI *R) ! "Syy" 

ASTRS (L, 3 ) = -FF*PR*COSD (THETAD) / (PI *R) ! "Szz" 

ASTRS ( L, 4 ) = 0 ! "Syz" 

ASTRS (L, 5) = 0 ! "Sxz" 

ASTRS (L, 6) = -FF*COSD (THETAD) **2*SIND (THETAD) / (PI*R) ! "Sxy" 

ASTRN (L, 1) = (ASTRS(L, 1) -PR* (ASTRS (L,2) +ASTRS (L,3) ) ) /EX ! "Exx 

ASTRN (L, 2) = (ASTRS(L,2) -PR* (ASTRS (L, 1) +ASTRS (L,3) ) ) /EX ! "Eyy 

ASTRN (L, 3) = (ASTRS(L, 3) -PR* (ASTRS (L, 1) +ASTRS (L,2) ) ) /EX ! "Ezz 

ASTRN (L, 4 ) = (1+PR) /EX*ASTRS (L, 4 ) ! "Eyz" 

ASTRN (L, 5) = (1+PR) /EX*ASTRS (L, 5) ! "Exz" 

ASTRN (L, 6) = ( 1+PR) /EX* ASTRS (L, 6 ) ! "Exy" 

1310 CONTINUE 


DO 1320 I = 1,20 
R = COORD (1,1) 

THETA = COORD (1,2) *PI/180 

AD I S P ( 1 , 1 ) = -FF* ( 1- PR**2 ) * (LOG (R) +SIN (THETA) **2/ (2 * ( 1 -PR) ) )/ 

+ ( PI *EX) ! "ul" 

AD I S P ( 1 , 2 ) = -FF* (1 + PR) *( (1-2*PR) * THETA- SIN (THETA) *COS (THETA) ) / 
+ ( 2 *PI *EX) 1 "u2 " 

AD I S P ( 1 , 3 ) = 0 ! "u3 " 

1320 CONTINUE 
1400 RETURN 
END 


SUBROUT INE INTEG ( KR , KS , KE , VEDI , VI I , VI2 ) 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

* The subroutine calculates the terms necessary for local interaction * 

* integral and equivalent domain integral of each element . * 

* * 

* Subroutine input: KR , KS , KE * 

* Subroutine output: VEDI, VII, VI2 * 

■k k 
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Definition of local variables: * 

ADUXl(i,j) Array for the derivatives of the j-th auxiliary displacement * 
component in calculating Til, w.r.t. the local xl coordinate at* 
the i-th integration point. j=l:ul; j=2:u2; j=3:u3. * 

ADUX113(i,j) Array for the derivatives of the j-th auxiliary * 

displacement component in calculating T13, w.r.t. the * 

local xl coordinate at the i-th integration point. j=l:ul;* 
j=2:u2; j=3:u3. * 

ASIGMA(i,j) 3x3 array for the auxiliary stress tensor of the element in * 
calculating Til. * 

ASIGMA13 ( i , j ) 3x3 array for the auxiliary stress tensor of the element * 

in calculating T13 . * 

Coefficients for each of the 6 stress-work density terms in the* 
calculation of the interaction integral. * 


* Cl ( i) 
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CEDI (i) 
DSHAP ( i , j ) 
DSRST ( i , j ) 

DSX ( i , j ) 

DUX1 ( i , j ) 

KE 

KR 

KS 

RX ( i , j ) 
SF(i) 

SIGMA (i, j) 
TERMED I (i, j 

TERMI2 (i, j) 

TERMII (i, j) 

WEDI (k) 

WI2 (k) 

WII (k) 

VEDI (i, j , k) 

VI2 (i , j , k) 

VII (i, j ,k) 


XRJ 

k k k k k k 


Coefficients for each of the 6 stress-work density terms in the* 
calculation of the equivalent domain integral. * 

Array for the derivative of the i-th shape function w.r.t. the * 
local coordinate j. j=l:xl; j=2:x2; j=3:x3. * 

Array for the derivatives of the s-function w.r.t. the natural * 
coordinates j at the i-th integration point. j=l:xi; j=2:eta; * 
j =3 : zeta. * 

Array for the derivatives of the s-function w.r.t. the local * 

coordinates j at the i-th integration point. j=l:xl; j=2:x2; * 

j=3:x3. * 

Array for the derivatives of the j-th displacement component * 
w.r.t. the local xl coordinate at the i-th integration point. * 
j=l:Ul; j=2:u2; j=3:u3. * 

The element number in the ring. * 

The ring number. * 

The element layer (segment ) number along the crack front. * 

Inverse of the Jacobian. i=1..3; j=1..3 * 

s-function at the i-th integration point, i=1..8. * 

3x3 array for the stress tensor of the element. * 

) Array for the terms in the expression of the equivalent * 

domain integral. i=l : 1st term; i=2 : 2nd term. * 

Array for the terms in the expression of the interaction * 

integral 1(2). i = l : 1st term; i=2 : 2nd term. * 

Array for the terms in the expression of the interaction * 

integral 1(1) . i=l : 1st term; i=2 : 2nd term. * 

Array for the stress- work density at the k-th integration point* 
for equivalent domain integral. * 

Array for the stress-work density at the k-th integration point* 
in the calculation of the interaction integral 1(2) . * 

Array for the stress-work density at the k-th integration point* 
in the calculation of the interaction integral 1(1) . * 

Array for the value of the equivalent domain integral of the * 
k-th element in the (i+l)-th ring of the j-th layer (segment) . * 

Array for the value of the interaction integral 1(2), for T13 * 

stresses, of the k-th element in the (i+l)-th ring of the j-th * 
layer (segment) . * 

Array for the value of the interaction integral 1(1), for Til * 
stresses, of the k-th element in the (i+l)-th ring of the j-th * 
layer (segment) . * 

Determinant of the Jacobian. * 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON PI , PL , HEAD 

COMMON /AUX/ ASTRS(8,6) ,ASTRN(8,6) ,ADISP(20,3) 

COMMON / AUX1 3 / ASTRS13 (8,6) , ASTRN13 (8,6) , ADISP13 (20,3) 

COMMON /FESOL/ DISP (2 1 , 3 ) , STRS ( 8 , 6 ) , STRN ( 8 , 6 ) 

COMMON /FEMOD/ MIE ( 50 , 24 ) , MNE ( 50 , 24 , 2 0) , COORD (2 1 , 3 ) , COORDL ( 2 0 , 3 ) 
COMMON /GIP/ RG ( 8 ) , SG ( 8 ) , TG ( 8 ) 

DIMENSION SF ( 8 ) , DSRST (8,3), DSX (8,3) , DUX1 (8,3), SIGMA (3,3) 
DIMENSION DSHAP (20,3) , RX (3,3) 

DIMENSION AS IGMA (3,3) ,ASIGMA13 (3,3) ,ADUX1 (8,3) ,ADUX113 (8,3) 
DIMENSION CEDI (6) , TERMEDI (8,2) ,VEDI (4,50,24) , WEDI (8) 

DIMENSION Cl (6) , TERMII (8, 2) , VII (4 , 50 , 24 ) , WII (8) 

DIMENSION TERMI2 (8,2) ,VI2 (4,50,24) , WI2 (8) 

DATA CEDI /0. 5, 0.5, 0.5, 1,1,1/ 

DATA Cl /1,1,1,2,2,2/ 


= 0 
= 0 


C. . . Initialization: 
VEDI (KR, KS , KE) 
VII ( KR , KS , KE ) 
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VI 2 ( KR , KS , KE ) = 0 
DO 1110 K = 1, 8 
WEDI (K) = 0 
WII(K) = 0 
WI2 (K) = 0 
TERMEDI ( K, 1 ) = 0 
TERMII(K,1) = 0 
TERMI2 (K, 1) = 0 
DO 1105 KI = 1,3 
DSX (K, KI) = 0 
DUX1 (K, KI) = 0 
ADUX1 (K, KI) = 0 
ADUX113 (K, KI) = 0 
1105 CONTINUE 
1110 CONTINUE 

C... DO-Loop #1199 - Sum over 8 integration points: 

DO 1199 LI = 1,8 ! LI : counter for the integration point 

CALL SHPF (DSHAP , XRJ, RX, LI , RG, SG,TG) 

C... Define s-function and its derivatives: 

SF (LI) = 0.5* (1+SG(LI) ) * (l-TG(LI) **2) ! s (xi , eta , zeta) 

DSRST ( LI , 1 ) = 0. ! ds/d (xi ) 

DSRST ( LI , 2 ) = 0.5* (l-TG(LI) **2) ! ds/d (eta) 

DSRST (LI, 3) = -1* (1+SG(LI) ) *TG(LI) ! ds/d(zeta) 

DO 1120 KI = 1,3 
DO 1120 KJ = 1,3 

DSX ( LI , KI ) = DSX (LI , KI ) +RX ( KI , KJ) *DSRST ( LI , KJ) 

1120 CONTINUE 

C ... Calculate the stress-work density W for the equivalent domain integral and 
C the interaction integral: 

DO 1130 K = 1,6 

WEDI (LI) = WEDI (LI) +CEDI (K) *STRS (LI,K) *STRN(LI,K) 

WII(LI) = WII(LI) +C1 (K) *STRS(LI,K) *ASTRN(LI,K) 

WI2 (LI) = WI2 (LI) +C1 (K) *STRS (LI, K) *ASTRN13 (LI,K) 

1130 CONTINUE 

C ... Calculate [du/dxl] and [du/dxl]_a terms: 

DO 1140 KI = 1,3 

DO 1135 KN = 1,20 ! KN : counter for the local nodes 

DUX1 (LI , KI) = DUX1 (LI , KI ) +DSHAP ( KN , 1 ) *DISP ( KN , KI ) 

ADUX1 ( LI , KI ) = ADUX1 (LI , KI) +DSHAP (KN, 1 ) *ADISP (KN, KI ) 

ADUX113 (LI, KI) = ADUX113 ( LI , KI) +DSHAP ( KN, 1) *ADISP13 ( KN , KI ) 

1135 CONTINUE 
1140 CONTINUE 

C ... Construct the stress tensors: 

DO 1150 KI = 1,3 
DO 1145 KJ = 1,3 
IF (KI .EQ. KJ) THEN 
S IGMA ( KI , KJ) = STRS (LI , KI) 

ASIGMA (KI , KJ) = ASTRS (LI , KI) 

ASIGMA13 (KI, KJ) = ASTRS13 (LI , KI ) 

ELSE 

SIGMA (KI , KJ) = STRS (LI, 9-KI-KJ) 

ASIGMA (KI , KJ) = ASTRS (LI, 9-KI-KJ) 

ASIGMA13 (KI , KJ) = ASTRS13 (LI , 9-KI-KJ) 

END IF 

1145 CONTINUE 
1150 CONTINUE 

C. .. Calculate the first terms of the equivalent domain integral and the 
C interaction integrals: 

C EDI- (du/dxl) *sigma* (dS/dx) <TERMEDI (LI , 1) > 


! K: counter for the integration points 


! KI : counter for the local coordinates xl to x3 
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C II - (du/dxl)_a*sigma* (dS/dx) + (du/dxl) *sigma_a* (dS/dx) <TERMII (LI , 1) > 
DO 1160 KI = 1,3 
DO 1160 KJ = 1,3 
TERMEDI (LI, 1) = TERMEDI (LI , 1) + 

+ DUX1 (LI , KI) *SIGMA (KI , KJ) *DSX(LI,KJ) 

TERMII ( LI , 1 ) = TERMII ( LI , 1) + (ADUX1 (LI , KI ) *SIGMA (KI , KJ) + 

+ DUX1 (LI , KI) *ASIGMA (KI , KJ) ) *DSX(LI,KJ) 

TERMI2 (LI , 1) = TERMI2 ( LI , 1) + (ADUX113 (LI , KI ) *SIGMA (KI , KJ) + 

+ DUX1 (LI , KI) *ASIGMA13 (KI,KJ) ) *DSX(LI,KJ) 

1160 CONTINUE 

C ... Calculate the second terms of the equivalent domain integral and the 
C interaction integrals - W*(dS/dxl) : 

TERMEDI (LI, 2) = WEDI (LI) *DSX (LI , 1 ) 

TERMII (LI , 2 ) = WII (LI) *DSX (LI , 1) 

TERMI2 (LI , 2 ) = WI2 (LI) *DSX (LI , 1) 

C ... Calculate the integrals for the element: 

VEDI (KR, KS , KE) = VEDI ( KR , KS , KE) + (TERMEDI (LI , 1 ) -TERMEDI (LI , 2 )) *XRJ 
VII (KR,KS,KE) = VII (KR,KS,KE) +(TERMII (LI, 1) -TERMII (LI, 2) ) *XRJ 
VI 2 (KR, KS , KE) = VI2 (KR, KS, KE) + (TERMI2 (LI , 1) -TERMI2 (LI , 2) ) *XRJ 

1199 CONTINUE 

1200 RETURN 
END 


SUBROUT INE SHPF (DSHP, XSJ, SX, L, RG, SG, TG) 

**************************************** 


* The subroutine forms the shape functions and their derivatives for the * 

* 20-node 3-D solid element. The orientation of the local nodes 1 to 20 is * 

* based on the ANSYS SOLID95 element type. * 

* Ref.: I.M. Smith & D.V. Griffiths, Programming the Finite Element Method, * 

* pp. 432-433. John Wiley & Sons, 1988. * 

* ★ 


* Subroutine input: L , RG , SG , TG * 

* --- Subroutine output: DSHP, XSJ, SX * 

* * 


* 

•k 

* 

* 

★ 

★ 

★ 

★ 

★ 

* 

* 


Definition of local variables: * 

DER(i,j) Array for the derivative of the i-th shape function w.r.t. the * 
natural coordinate j. j=l:xi; j=2:eta; j=3:zeta. * 

FUN(i) Array for the shape function of the i-th node w.r.t. the * 

natural coordinates. * 

IR ( i) , IS ( i) , IT ( i) Arrays for the natural coodinates of the i-th node. * 

L Counter for the integration points. * 

SX(i,j) The inverse of the Jacobian. * 

XS(i,j) The Jacobian. * 

XSJ The determinant of the Jacobian. * 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /FEMOD/ MIE ( 50 , 24 ) , MNE ( 50 , 24 , 2 0) , COORD (2 1 , 3 ) , COORDL ( 2 0 , 3 ) 
DIMENSION IS(20) ,IT(20) ,IR(20) ,SG(8) ,TG(8) ,RG(8) 

DIMENSION FUN (20) ,DER(20,3) ,XS(3,3) ,SX(3,3) ,DSHP(20,3) 

DATA IR /-1,1, 1,-1, -1,1, 1,-1, 0,1, 0,-1, 0,1, 0,-1, -1,1, 1,-1/ 

DATA IS /-l, -1,1, 1,-1, -1,1, 1,-1, 0,1, 0,-1, 0,1, 0,-1, -1,1,1/ 

DATA IT /-l, -1,-1, -1,1, 1,1, 1,-1, -1,-1, -1,1, 1,1, 1,0, 0,0,0/ 


C... Define shape functions and their derivatives for each node: 

DO 1205 I = 1,20 ! I: counter for the local node numbers 

R = RG (L) *IR ( I) 

S = SG (L) *IS(I) 

T = TG (L) *IT(I) 

IF ( IR ( I ) .EQ. 0) THEN ! local nodes 9,11,13,15 

FUN ( I ) = 0.25* (l-RG(L) **2) * (1+S) * (1+T) 
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DER(I,1) = -0.5*RG(L) * (1+S) * (1+T) 

DER(I,2) = 0.25*IS(I) * (l-RG(L) **2) * (1+T) 

DER(I,3) = 0.25*IT(I) * (l-RG(L) **2) * (1+S) 

ELSEIF (IS(I) .EQ. 0) THEN ! local nodes 10,12,14,16 

FUN ( I ) = 0.25* (1+R) * (l-SG(L) **2) * (1+T) 

DER (1,1) = 0.25*IR(I) * (l-SG(L) **2) * (1+T) 

DER (1,2) = -0.5*SG(L) * (1+R) * (1+T) 

DER (1,3) = 0.25*IT(I) * (1+R) * (l-SG(L) **2) 

ELSEIF ( IT ( I) .EQ. 0) THEN ! local nodes 17,18,19,20 

FUN ( I ) = 0.25* (1+R) * (1+S) * (l-TG(L) **2) 

DER (1,1) = 0.25*IR(I) * (1+S) * (l-TG(L) **2) 

DER (1,2) = 0.25*IS(I) * (1+R) * (l-TG(L) **2) 

DER (1,3) = -0.5*TG(L) * (1+R) * (1+S) 

ELSE ! local nodes 1,2,3,4,5,6,7,8 

FUN ( I ) = 0.125* (1+R) * (1+S) * (1+T) * (R+S+T-2) 

DER (1,1) = 0. 125*IR(I) * (1+S) * (1+T) * (2*R+S+T-1) 

DER (1,2) = 0. 125*IS (I) * (1+R) * (1+T) * (R+2*S+T-1) 

DER (1,3) = 0. 125*IT(I) * (1+R) * (1+S) * (R+S+2*T-1) 

END IF 

1205 CONTINUE 

C ... Construct the Jacobian, its determinant, and the inverse of the Jacobian: 

DO 1210 I = 1,3 
DO 1210 J = 1,3 
XS(I,J) = 0 
DO 1210 K = 1,20 

XS(I,J) = XS (I , J) +COORDL (K, J) *DER (K, I) ! Jacobian 

1210 CONTINUE 

CALL MINV ( SX , XS J , XS ) 

C...Form derivatives of the shape functions in global coordinates. 

DO 1230 I = 1,20 
DO 1230 J = 1,3 
DSHP ( I , J) = 0 
DO 1230 K = 1,3 

DSHP ( I , J) = DSHP (I , J) +SX ( J, K) *DER ( I , K) 

1230 CONTINUE 

1300 RETURN 
END 

SUBROUTINE MINV (AINV, DET, A) 

* The subroutine calculates the determinant of a 3x3 matrix, and forms its* 

* inverse matrix. A standard Gauss-Jordan elimination algorithm is used. * 

* Ref.: G.J. Borse, FORTRAN 77 and Numerical Methods for Engineers, pp.429- * 

* 432. PWS Publishers, 1985. * 

* * 

* Subroutine input: A * 

* Subroutine output: AINV, DET * 

■k k 

* Definition of local variables: * 

* A The input matrix. * 

* AINV The inverse of A. * 

* DET The determinant of A. * 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION A (3, 3) , AINV (3, 3) ,B(3,3) 

DO 1 I = 1,3 
DO 1 J = 1 , 3 
B (I, J) = A ( I , J) 

IF (I .EQ. J) THEN 
AINV ( I , J) = 1 
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ELSE 

AINV(I,J) = 0 
ENDIF 

1 CONTINUE 

DO 10 IPASS =1,3 ! IPASS: counter for the current pivot row 

C. . .For each pass, find the maximum element in the pivot row: 

IMX = IPASS 

DO 2 IROW = IPASS, 3 

IF (ABS (B( IROW, IPASS) ) .GT. ABS (B (IMX, IPASS) ) ) THEN 
IMX = IROW 
ENDIF 

2 CONTINUE 

C. . .Interchange the elements of row IPASS and row IMX in both B and AINV : 
IF (IMX .NE. IPASS) THEN 
DO 4 ICOL = 1,3 
TEMP = AINV (IPASS, ICOL) 

AINV (IPASS, ICOL) = AINV (IMX, ICOL) 

AINV (IMX, ICOL) = TEMP 
IF (ICOL .GE. IPASS) THEN 
TEMP = B (IPASS, ICOL) 

B (IPASS, ICOL) = B (IMX, ICOL) 

B (IMX, ICOL) = TEMP 
ENDIF 

4 CONTINUE 

ENDIF 

PIVOT = B (IPASS, IPASS) ! the current pivot 

C. . .Normalize the pivot row by dividing across by the current pivot: 

DO 6 ICOL = 1,3 

AINV (IPASS, ICOL) = AINV (IPASS, ICOL) /PIVOT 
IF (ICOL .GE. IPASS) THEN 
B( IPASS, ICOL) = B( IPASS, ICOL) /PIVOT 
ENDIF 

6 CONTINUE 

C. . .Replace each row by the row plus a multiple of the pivot row with the 
C factor chosen so that the element of [B] in the pivot column is 0: 

DO 8 IROW = 1,3 
IF (IROW .NE. IPASS) THEN 

FACTOR = B( IROW, IPASS) ! set the factor for this row 

ENDIF 

DO 7 ICOL = 1,3 
IF (IROW .NE. IPASS) THEN 

AINV( IROW, ICOL) = AINV(IROW, ICOL) -FACTOR *AINV ( IPASS, ICOL) 

B (IROW, ICOL) = B (IROW, ICOL) -FACTOR*B( IPASS, ICOL) 

ENDIF 

7 CONTINUE 

8 CONTINUE 
10 CONTINUE 

DET = A ( 1 , 1) *A (2,2) *A (3,3) -A(l,l) *A (2,3)*A(3,2) + 

+A (1,2) *A (2,3) *A (3,1) - A ( 1 , 2 ) *A (2,1) *A (3,3)+ 

+A (1,3) *A (2,1) *A (3,2) -A(l,3) *A (2,2) *A (3,1) 

RETURN 

END 


SUBROUT INE S TRN3 3 ( E 3 3 , KE ) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* The subroutine calculates the average Ezz strain of a wedge-shaped * 

* element from the finite element result. * 

■k k 


Subroutine input: KE, data file 8 
Subroutine output: E33 
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* Definition of local variables: * 

* KE The element number in the layer (segment ) . * 

* E33(i) Array for the average Ezz strain of the i-th element. * 

* STRNZZ (i,j) Array for the FE result of strains within a wedge-shaped * 

* element, which is attached on the crack front. The i-th row * 

* stores strain components at the i-th Gaussian integration * 

* point, i=1..8. Strain components from 1st to 6th column are * 

* [Exx, Eyy, Ezz , Exy , Eyz , Exz] , respectively. * 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION STRNZZ (8, 6) ,E33 (24) 

C...Read strain data within a wedge-shaped element (the element attached on the 
C crack front) : 

READ (8,1405) (STRNZZ ( 1 , J) , J=1 , 6 ) 

1405 FORMAT (////9X, 6 (E12 . 5) ) ! ANSYS 5.5 format 

1406 FORMAT ( 9X, 6 (E12 . 5 ) ) 

DO 1410 I = 2,8 

READ (8,1406) (STRNZZ ( I , J) , J=1 , 6) 

1410 CONTINUE 

C ... Calculate average strain_33 at the mid-side node on crack front: 

E33 (KE) = ( STRNZZ (4,3) + STRNZZ ( 8 , 3 ) ) /2 
1500 RETURN 
END 



